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EXECUTIVE SUMMARY 


This report supplements a previous report of the same title submitted in June, 1992. It 
summarizes additional analytical techniques which have been developed for predicting the 
response of linear and nonlinear structures to noise excitations generated by large propulsion 
power plants. 

The report is divided into nine chapters. The first two deal with incomplete knowledge 
of boundary conditions of engineering structures. The incomplete knowledge is characterized by 
a convex set, and its diagnosis is formulated as a multi-hypothesis discrete decision-making 
algorithm, with attendant criteria of adaptive termination. 

Chapter 1 deals with rectangular plates. The boundary conditions are represented by 
transverse and rotational springs with uncertain spring constants, modeling possible damages 
along the boundary. Bolotin’s dynamic edge effect method is applied to determine the natural 
frequencies. Chapter 2 is devoted to the identification of end conditions for beams, and their 
natural frequencies, using a finite element formulation. 

Chapter 3 deals with free and forced vibrations of periodic multi-span beams. The 
concept of wave propagation in periodic structures of Brillouin is utilized to investigate the wave 
motion at the periodic supports. Chapter 4 discusses the vibrations of a two-dimensional grillage. 
Each periodic support of such a grillage is constrained by a rigid transverse support and two 
elastic rotational springs in two orthogonal directions. The elastic springs on each row are 
identical; thus, the grillage forms a two-dimensional periodic pattern. The four boundary edges 
of the grillage are assumed to be either simply supported or clamped. Again the wave 
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Propagation concept is used for the analysis. 
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Diagnosis of Local Modifications in 

the Boundary Conditions of a 
Rectangular Plate 
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DIAGNOSIS OF LOCAL MODIFICATIONS IN THE 
BOUNDARY CONDITIONS OF A RECTANGULAR PLATE 


Isaac Elishakoff and Jianjie Fang 

Center for Applied Stochastics Research 
and Department of Mechanical Engineering 
Florida Atlantic University 
Boca Raton, FL 33431-0991, USA 

1. INTRODUCTION 

There is a vast literature on the free vibrations, as well as forced vibrations, either due 
to deterministic or random excitations, of thin rectangular plates. To the best of our knowledge, 
in all these studies the boundary conditions along the edges of the plates are assumed as known. 
Hence the investigators consider all possible combinations of boundary conditions along edges, 
namely those of simply supported (SS), clamped (C), free (F), or, more generally, elastically 
supported (ES) conditions. The natural question arises: How to identify the "true" boundary 
conditions? In addition, the integrity of plate-like structures is often determined by the condition 
at the junction of the plate with the remainder of the structure. Diagnosis of the integrity of the 
junction can be performed by estimating the existing boundary conditions. 

In this study the boundary condition diagnosis is formulated as a discrete multi-hypothesis 
decision problem. This facilitates diagnostic determination of global features of the boundary 
conditions, without the necessity of detailed full reconstruction of the boundary condition profile 
along the entire perimeter of the plate. In particular, we concentrate on diagnosis of localized 
alteration of the boundary conditions due to some adverse effects, including the high level 
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specified boundary conditions, which can be performed analytically or numerically. This is 
necessary for implementing the minimum-distance decision algorithm upon which the diagnosis 
is based. In the present study we employ the generalized Bolotin’s dynamic edge effect method 
to determine the approximate natural frequencies and normal modes of elastically supported 
isotropic, uniform rectangular plates. The entire procedure of diagnosis of local modifications 
is implemented with numerically simulated diagnosis. 

2. CONVEX MODELS OF UNCERTAINTY 

As indicated above, diagnosis of the boundary condition means, in this study, 
identification of the affected region on the boundary and estimation of limits on the magnitude 
of the stiffness change. It is assumed that only the torsional stiffness is modified by boundary 
condition failure. Convex models put forward by Ben-Haim and Elishakoff [1] will be used to 
specify the degree of precision required in the diagnosis of the boundary torsional stiffness 

profile. 

A convex model is a set of functions. Each element of the set represents a possible 
realization of an uncertain quantity of interest. In our case, these functions are the torsional 
stiffness profiles k(s), where s is the position around the periphery of the plates. The diagnosis 
will be considered satisfactory when the uncertainty in the boundary condition has been reduced 
to an acceptable preselected level. Convex models will represent the acceptable uncertainty in 
the boundary condition. This means that diagnosis of the boundary condition is in fact no more 
than identification of the convex model to which the actual boundary stiffness profile belongs. 
It should be stressed again that the actual stiffness profile will not be identified; only the convex 
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model to which it belongs is sought by the diagnosis. 

While several different convex models are available (see Ben-Haim and Elishakoff, [1]), 
the following simple model is particularly suitable for representing uncertain localized failures. 
Let S represent the boundary of the plate and let /„ , n = denote subsets of the boundary 

5. The sets /„ may be simple connected intervals, or unions of such intervals. Let k(s) be the 

nominal torsional stiffness profile, s £ S. We define convex models which are sets of stiffness 
profiles which deviate from the nominal profile only in a particular interval Specifically, the 
convex model is the following set of hypothesized stiffness profiles [2]: 

C„ = : Ks) - m, *£/.; 1*0 - 1*6. (1) 

In other words, is the set of torsional stiffness profiles k(s) which are equal to the 
nomi nal function k(s) outside the region /„, but which deviate by as much as ±6 from the 
constant value K m throughout the region /„. One recognizes that C m is a localized uniform bound 
convex model 

3. MULTI-HYPOTHESIS DECISION 

The reference values K m in the convex models in Eq. (1) can assume any of the M 
different values K v Likewise, failure (i.e. "weakening" of the boundary conditions 

intended in the original design) can occur in any of N different boundary intervals, l v ..., I N - 
Thus diagnosis of the boundary condition involves deciding which of the MN convex models 
contains the true boundary stiffness profile. This decision is based on a multi-hypothesis 
formulation. 


1-4 


In its simplest form the multi-hypothesis decision algorithm requires the choice of one 
representative or hypothesized stiffness profile from each convex model. It should be noted that 
in its more general form, more than one hypothesized stiffness profile is chosen for each convex 
model. Let hj^s) denote the hypothesized stiffness profile from convex model C m . That is, 
h m»( s ) e c mn- Let indicate a vector of / natural frequencies of a plate 

whose boundary stiffness profile is h^s). Furthermore, let Q = (Q Jt ...£2/) represent the vector 
of measured natural frequencies of the same modes. Let ||x| denote a norm of the vector x. 

The distance from the measured natural frequency vector Q to the anticipated natural 
frequency vector based on the mnth hypothesis is : 

H m = | Q - 

It is reasonable to conclude that the true stiffness profile is likely to belong to the convex 
model whose anticipated spectrum of natural frequencies is closest to the measured spectrum. 
Let the pair of indices m 0 and n 0 satisfy: 


H 


<v. 


= min H_ 

m,n 


(3) 


The multi-hypothesis decision is that the true stiffness profile belongs to the convex model 


C . 




4. DIAGNOSIS PROCEDURE 

The decision algorithm implied by Eq. (3) will always reach a decision, regardless 
whether or not the least distance, H is small. In other words, the multi-hypothesis will choose 
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the most likely one from among the available options, but none may be particularly convincing, 
due to unsatisfactory choice in the preliminary analysis. 

Based on the assumption that all needed low-order frequencies of the plate are measured, 
the diagnosis procedure which we follow, in broad outline, is as follows: 

1. Construct the convex model of the unknown boundary conditions by taking advantage 
of the available prior knowledge of the range of values which can reasonably be expected to 
occur. If such a knowledge is unavailable, 

2. Calculate all the natural frequencies of interests for each hypothesized stiffness profile 
hmn( s ) by employing the generalized Bolotin’s dynamic edge effect method, or any other 
approximate, analytical, or purely numerical technique. 

3. Compute the distances between the measured frequencies and the ones corresponding 
to hypothesized stiffnesses according to the general definition given in Eq. (2). 

4. Decide which convex model is the most likely one, on the basis of minimum distance 
algorithm relation (3). If is larger than a specified threshold value 0 then decision is deferred, 
and a new choice of the mostly convex model is made. The diagnosis process is terminated if 

H * 0 

mn 

As is seen, one of the cornerstone of the method is the ability to determine the natural 
frequencies of plates under arbitrary boundary conditions. This problem is addressed in the 
following section. 


5. CALCULATION OF EIGENFREQUENCY SPECTRA 
There exists an ample literature on vibration of elastic plates, as described in definitive 
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reviews by Leissa [3-5]. For plates with two opposite edges simply supported but with arbitrary 
boundary conditions along the remaining edges, exact solution is available. This class of 
vibration problems is referred to as Livy problems. For non-Livy problems exact solutions are 
presently unavailable and approximate techniques should be resorted to. Energy based 
techniques, like methods of Rayleigh, Rayleigh-Ritz, Galerkin, or the finite element methods are 
accurate and relatively incumbersome for the determination of the lower end of frequency 
spectrum. However, for higher frequencies use of energy methods becomes computationally 
expensive and cumbersome. 

For higher frequencies, Bolotin [6] proposed so called dynamic edge effect method. It 
is based on the physical observation, that for high frequencies the mode shapes of plates of 
different boundary conditions have, a somewhat similar behavior except in the zones near edges, 
where the effect of the boundary conditions is prominent. Bolotin therefore suggested to 
construct solutions emanating from each edge. In the interior region he "tailored" these solutions 
to yield two simultaneous characteristic equations. This method found ample popularity both in 
the East and the West (for appropriate bibliography one may consult review article by Elishakoff 
[7] ). It turned out, however, that for certain instances it was impossible to construct solutions, 
having sufficient number of decaying terms, to converge to an interior solution. This 
phenomenon was called by Bolotin degeneracy of the dynamic edge effect method. The "cure" 
for such a situation was provided by Vijayakumar [8] and Elishakoff [9]. In particular in Ref. 
9 it was suggested to solve two auxiliary Levy-type problems in conjunction with the relationship 
between the natural frequency and the modal numbers. In this formulation Bolotin’s idea of 
necessity of decaying solutions towards the interior region was abandoned, and solution was 
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required to be valid in hypothesized "strips", somewhat similar to the original idea of Germaine 
(modified by Lagrange) to derive the very equations for the plate. Here we follow the 
generalization of Bolotin’s dynamic edge effect method for the elastically supported plates, but 
in the framework of Ref. 9. 

Differential equation of motion of Kirchoff-Love plate reads 

DV 2 Vhv + ph— = 0 ( 5 ) 

dt 2 

where D is the flexural stiffiiess, p-material density, A-thickness, w(x,y,t) - displacement 
V 2 = d 2 /dx 2 + d 2 ldy 2 Laplace operator, x and y space coordinates, t - time. The boundary 
conditions for rotationally restrained rectangular plate are 


dw 

= a ‘17 


x = 0 


~ 'i'Mx = a : 


dw 

dx 


x - a 


( 6 ) 


(7) 


sm = 

1 y Kl dy 


y =0 


( 8 ) 


w = 0 on S 

where a and b are side lengths of the plate, 


(9) 

( 10 ) 
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( 11 ) 


II 

b 

dhv 

d 2 w\ 

+ v 

X j 

[dx 2 

By 2 ) 


( 


M - D 

d 2 * 

d 2 w 

+ V — 

y 


dx 2 j 


( 12 ) 


are the bending moments per unit length. Coefficients a 2 and a 2 are stiffnesses of rotational 
springs in the x direction; and p 2 are stiffness of rotational springs in the y direction. In 
addition, artificial parameters y ly and 6 2 are introduced. They take values either zero or 
unity, in order to accommodate various possible ideal boundary conditions. For example y 2 = 0, 
a 2 * 0 corresponds to the edge x - 0 being clamped; combination 0, a,= 0 corresponds to 
simply supported edge. 

For free vibrations we seek solution of the Eq. (5) in the form 

Mx,y,t) = W(x,y)e “* ^ 

where W(x,y) is the vibration mode and w is the sought natural frequency. Eq. (5) becomes 


Z>V 2 V 2 W - p h<a 2 W = 0 


(14) 


Consider first the plate which is simply supported all round S, with boundary conditions 


W = 0 , V 2 W = 0 
The mode shape 


(15) 


W - sin^Lsin^L , (u, v= 1,2,— ) 

“ v a b 


(16) 


satisfies both the governing differential equation (14) and the boundary conditions (15). The 
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serial numbers u and v are positive integers, representing, respectively, number of half waves in 
x - and - y - directions. The substitution of Eq. (16) into Eq. (14) yields the expression of the 
natural frequency 


2 

CO = 
uv 


D 
P h 




We introduce the nondimensional frequency parameter X^, defined by 


(17) 


= u) w a 2 s[ph}D = Jt 2 [u 2 * ( a/b) 2 v 2 ] ( 18 ) 

Consider now the plate with original boundary conditions specified in Eqs (6-10). 
Following the spirit of the Bolotin’s dynamic edge effect method [6] as well as its generalization 
by Elishakoff [9], we represent the relationship between the natural frequency squared and modal 
numbers in a manner similar to the Eq. (17) 


(19) 

where under new circumstances p and q, are unknown positive decimal numbers rather than 
integers. Similar to Eq. (18), the corresponding nondimensional frequency X is 


co 2 = 


ph 


(pn\ 

2 

+ 

' qn\ 

2 



\ b ) 



X^ = a)a 2 \jph/D = n 2 \p 2 +(a/b) 2 q 2 ] 

Eq. (14) becomes 


( 20 ) 
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V 2 V 2 W - 


' pn 


f qn'' 

It 


W = 0 


( 21 ) 


The solution of this equation can be represented as composed of solutions W 1 and W 2 of 
following problems 



( 22 ) 


(23) 


These will be determined, as suggested by Elishakoff [9], from two auxiliary problems 
discussed below. 


FIRST AUXILIARY PROBLEM 

Consider the first auxiliary problem, referred to as X-Problem. We seek for the solution 
of Eq. (21) in the form 

W(x,y) = X(x) sin^L (24) 

0 

with boundary conditions, (6), and (7) and (10). Instead of Equations (22) and (23) we arrive 
at 



( 25 ) 


1-11 


dx : 


’pn) 

a 


qn\ 2 


+ 2Z1 
b 


) J 


1 * 2=0 


(26) 


The solution for the function X(x) 


X(x) = X^x) + X 2 (x) 


(27) 


reads 


X(x) = AjSin£^f. + A 2 cosE^L + A 3 sinh.!ZH + A 4 cosh 


nvc 


(28) 


where the first two terms are associated with Eq. (25) and the last two terms are associated with 
Eq. (26). In Eq. 26 the following notation has been used 


r = p 


/ — \ 2 / _ \2 

1 +2fl] w 

[P [b 


and Aj arc constants of integration. 

Boundary conditions (6), (7) and (10) lead following conditions for X(x) 
X = 0, at x = 0, x - a 


(29) 


(30) 


„d 2 X dX 
Y £> = a. — at 

1 dx 2 l dx 


x = 0 


(31) 


„d 2 X dX . v „ 

l^-TT = a 2~r~ at x = a 


(32) 


<£c 2 z d!x 

In writing Eqs (31) and (32) we have used the boundary conditions (6) and (7). 
Satisfaction of boundary conditions leads to: 
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(33) 


A 2 + A 4 = 0 

A l S l + Vl + A A + = 0 

A l <j» 1 + Ayjm + A y R^ l - A^jm = 0 

A i(+2 c i ~ " A 2 (<*»/ 1 + YaP^i) 

+ A 3 (<^ 2 i?C 1 + y^ 2 P kS i) + + r 2 P^i C i) = 0 

where 


(34) 

(35) 


(36) 


♦ i = 


a i a 


'I- — 


R - — 


D-D p 

= sin(/m) , c L = cos (pn) , S l = sinh(nt) , C x = cosh(m) 


(37) 


<^ L and <^ 2 are nondimensional rotational spring stiffnesses. 

Requiring that A ; ’s do not vanish simultaneously results in a frequency determinant of the 
X-problem 


0 


4>t 

4»2 c r YaP* 5 ! 


1 

C 1 

Yjpir 

_ (4 > 2 s i + YaP 7lc i) 


0 

Sx 


1 

-R 2 yj>Tt 


= 0 < 38 ) 


$ 2 RS l +R 1 pT:y 2 C l 


This equation contains two unknowns p and q. We will obtain a needed companion equation for 
determining these two unknowns by solving an analogous problem in the y direction. 
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SECOND AUXILIARY PROBLEM 


Consider now the second auxiliary problem referred to as an Y-Problem. We seek for the 
solution of Eq. (21) in the form 


W{x,y) = Y(y) sin£^ 
a 


(39) 


with boundary conditions (8) and (9). Instead of Equations (22) (23) and (10) we have 


£21 ♦/£ 

dy 2 \ b 


\ 2 


y i-o 


(40) 


d% 

dy 2 


(qn\ +2 pn 


I J 


y 2 =0 


(41) 


The solution for the function Y(y) 


Y(y) = Y t (y) * Y 2 (y) 


(42) 


is given by 


Y(y) = B^sin^L 
b 


+ B 7 cos + 5 3 sinh— + B 4 cosh^2i 
2 b 3 b 4 b 


(43) 


where 


z = q 



(44) 


Note that z is obtained from the expression (29) for r, by formal replacement of p by q, 
of q by p, of a by b, and of b by a. 

The boundary conditions (8), (9), (10) read 
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Y = 0, 


y = b 


( 45 ) 


at y = 0 , 


bp£L at y = 0 

1 dy ! 'dy 


-a p£L at y-b 

^ dy* z dy 


Satisfaction of the boundary conditions leads to following set of equations 

b 2 +b 4 =0 

B 1 S 2 + B 2 c 2 + B A + B A C 2 * 0 

B i qj 1 + Bjb^qii + - B£ 2 b Y qn - 0 

B^2 C 2 ~ b 2 ( P tS 2) ~ B 2^2 S 2 + 

+ B^ 2 ZC 2 + 6 ^cpiSJ * s 2 + Z 2 qjtb 2 Cj) = 0 

where 

*•■£*■ z = | 

s 2 = sin(^Jt) , c 2 = cos (qn), S 2 - sinh(zn), C 2 * cosh(zji) 

Moreover, and i|) 2 are nondimensional rotational constants. 

Requirement of nontriviality of the constants of integration yields the 
determinant 


( 46 ) 

( 47 ) 

( 48 ) 

( 49 ) 

( 50 ) 

( 51 ) 

( 52 ) 

following 
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0 


1 


0 


1 


y 2 c 2 -b 2 qns 2 -(qf 2 +b 2 qncj 


*2 C 2 
Z ^ 1 -Z 2 b y qK 

■tyjZC 2 +b 7 Z 2 qnS 2 ty^S 2 +Z 2 qnb 2 C 2 


(53) 


In contrast to classical free vibration problems [3], here two equations in terms of 
determinants, namely Eq. (38) and Eq. (53), should be solved simultaneously to determine the 
decimal modal numbers p and q. Their substitution into Eq. (19) yields the desired natural 
frequencies. It should be borne in mind that in contrast to energy based methods the 
determination of higher natural frequencies does not constitute a more difficult task than that of 
the low natural frequencies. In fact, the numerical effort for numerical evaluation of any natural 
frequency is essentially identical. The additional advantage. of the present generalization of the 
Bolotin’s method consists in the possibility of evaluation of any preselected natural frequency. 
This is important in particular due to the diagnostic method adopted in section 3. The validity 
of the present generalization of Bolotin’s dynamic edge method is verified by comparing the 
present results with the results due to Leissa [4]. The comparisons are tabulated in the Appendix. 


6. NUMERICAL EXAMPLE 

1) Example 1 

For simplicity, we first concentrate on the boundary condition diagnosis of a square plate 
( a=b ) with three edges clamped and the fourth edge elastically supported (C-C-C-ES) (see Fig. 
1). The uncertain rotational stiffness of the latter edge of the plate, (J 2 , is represented by the 
convex model, given in Eq. (1). We need some partial information on the plate boundary 
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condition, i. e. on p 2 , to proceed with the diagnosis. Since no experimental results were available 
to present investigators, the simulated situation was constructed. 

To simulate possession of some information about the to-be-diagnosed plate, we proceed 
as follows. We first investigate the extreme case with all four edges being clamped. The natural 
frequency parameter of the fundamental frequency by the present approximate analytical method 
is \j=35.113. The plate under consideration is an aluminum plate (p -268.7 kg/m*) with both 
side lengths fixed at a=b=10 m, and thickness at h=0.1 m. We first determine, by numerical 

evaluation of the generalized Bolotin’s method, the value of the rotational stiffness such that 

the first mode frequency will be some fraction (say, 97%) of the corresponding first mode 
frequency of the all-round clamped plate, namely 'k I =34.060. Numerical search yields the value 

of the nondimensional rotational stiffness fJ 2 = l-522xl0 7 N. We note that in this case the 

numerical solution of Eq. (36) and Eq.(51) gives the values p=1.3417 and q=1.2848. For values 
of f} 2 greater than p 2 the natural frequencies will be identified with those of essentially C-C-C-C 

plate. One would expect that beyond this value of the rotational spring coefficient identification 
of the "true" boundary condition would be difficult since all the plates with various (5> p 2 

actually are the same. 

As it was mentioned above, performing actual experimental measurements was beyond 
the scope of the present study. In order to simulate the actual situation, we suppose that we may 
get, in principle, all needed measured frequencies of the plate whose rotational stiffness 
coefficient f} 2 along the fourth edge should be identified. Let us visualize that the set of 
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measured nondimensional natural frequencies is given in Table 1. The notation where i is 
the sequence number of the frequency, is introduced in addition to the two index notation of 
frequency X n , because repeated frequencies are considered as a single frequency in this study. 
We note in passing that the frequency parameters in Table 1 were obtained by setting fi 2 =10 6 N. 
We visualize that this actual value of the rotational stiffness is unknown to us in our hypothetical 
experiment and we should identify the convex model which represents it in a best way. 


Table 1 The first six nondimensional frequency parameters of the to-be-diagnosed plate, 
simulating the results of measurements 


V) 

31.942 

*t2) 

64.119 

^3) 

71.088 


101.052 

\s) 

117.479 


130.355 


Since within the present procedure we assume that the alteration of the boundary 
conditions occurs along the entire length of the side, we set that N=4 in Eq. (1) since the plate 
has four sides. Accordingly the subsets of the boundary S, /„ n=l, ••• , 4, are cq, (3 ; , f$ 2 . Here 

the nominal values of torsional stiffness are denoted by c^, P 2 . As mentioned before, 

we assume three edges of the square plate to be clamped and the remaining edge (edge 4 in this 
case) to be elastically supported. Therefore, the nominal values of torsional stiffness 
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corresponding to the three clamped edges, denoted by a t , P p a 2 , are infinite, that is 


a 2 = = Oj -* oo . The fact that the nominal value of the elastically supported edge, p 2 , is 


unknown and is to be identified suggests that it may need several estimations to reasonably 
represent the convex model for the stiffness P 2 . We denote these reference values P 2 where 

m indicates the number of estimated nominal values. Then the general convex model set in 
Eq.(l) becomes 

( 54 ) 

C„ 4 - (tvM-1.2: a,*o, f = 1.2; P, -P, IP, - Pz„l * 


or, simply 


q. * IPr IP, - p,.i * »„> 


( 55 ) 


since only p 2 , the stiffness of the elastically supported edge, is considered to be unknown. 

The natural question arises: Is the convex modeling of the uncertain elastic supports 
along the fourth edge able to "uncover” the best convex model closest to this simulated situation? 

To answer the above question it appears to be instructive to consider another question: 
"How to reasonably construct the convex model?" Conceivably, if the convex model of stiffness 
P 2 is constructed that can include the true value of fi 2 , the answer is "yes". Otherwise it makes 
no sense to try to single out the "best" stiffness from the wrong ones. It is quite natural to state 
that the estimated stiffness is likely to be far from the true one if the true stiffness can not be 
exactly or approximately represented by the elements of the hypothesized convex model. This 
happens since even the closest stiffness element in the convex set can be quite different from the 
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true one. Thus it is necessary to develop guidelines as how to choose the proper convex models, 
namely those which will necessarily contain possible realizations of an uncertain quantity of 
interest, and hence able to represent the uncertainty of the boundary conditions. 

To illustrate the procedure to construct the convex model, we will investigate a quite 
general example in this section. 

First we hypothesize the nominal value of stiffness P 2 , say P 2.1 =0-65xl0 6 AT, and also set 

6j=0.15x10 6 N. The according convex model is (0.5xl0 6 , 0.8xl0 6 )N. As we will see later, the 
diagnosed stiffness will be identified at the upper bound of the interval (see Table 2). This 
implies that either this upper limit is a true value of p 2 , or it is closest to the true p 2 . Since it 
is improbable to "hit" the true value at the first trial, we conclude that the present convex model 
is not comprehensive and needs to be improved. The fact that the diagnosed value of p 2 lies at 

the upper bound of the interval suggests that a larger reference value p^ x should be chosen. 

Next we select p^ 2 =2.5x10 6 N and set b 2 =Q.5xlQ 6 N. Then the corresponding interval becomes 

(2xl0 6 , 3x10 6 )N. As we can see from the results, the diagnosed stiffness p 2 in this case lies at 
the lower bound of the interval (see Table 3). Again we conclude that either this lower limit is 
the true value of p 2 or it is closest to p 2 within this "wrong" interval. This leads to creation of 
the final convex model (0.8xl0 6 , 2xl0 6 )Af (see Fig. 2), and we can be confident that it is able to 
include the true value of stiffness p 2 . In practical situations, we may need several estimated 
models to construct the "final" convex model. Obviously the above situation represents the 
typical case that may occur. The procedure, in detail, is as follows: 

1. Measure all the "necessary" natural frequencies of the actual plate (In this study, we assume 


1-20 



that they are all provided and listed in Table 1.) 

2. Construct a convex set of boundary conditions. In this example it is (0.5xl0 6 , 0.8xl0 6 )N. The 
lower bound of boundary stiffness interval is denoted by p 2 , and the upper one by p 2 . 

3. Calculate the natural frequencies of the plate with the boundary stiffness P 2 equal to the mid- 
point value of the interval, p 2 = 0.65xl0 6 N in this example. 

4. Compute the distances H between the measured and calculated natural frequencies with p 2 = p 2 

= 0.65xl0?N. It is well known that boundary conditions have a more pronounced effect in the 
range of the low frequencies than in the high frequency region. This means that we should put 
more "weight" when defining the norm on the low frequencies. Thus the distance between the 
natural frequencies in Eq. (2) can be defined specifically as: 

//„ - II a -W>| (56) 

»« 1 \ 

where A is the vector of nondimensionalized measured natural frequency Q, and lS” m) is the 
calculated one for the hypothesized stiffness. 

5. Compute the distance H at values p 2 = p 2 + A0, denoted by H($\ +A0) following the same 

way as step 4. We set Ap at l.OxlO 3 N. The derivative of the distance H with respect to p 2 can 
be approximately estimated by 

dWSb m #(p;+Ap)-ff(pl) (5?) 

d& " H(jfy 
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6. Decide search direction of p 2 by the following manner: 


a) If H( P 2 + pj ) < 0, the stiffness p 2 should be increased. 

b) If H( P2 + Pj ) > 0, the stiffness p 2 should be decreased. 

7. Choose the new search interval. In this study, it is always one half of the last search interval, 
and is created in such a way that the actual p 2 is always in the new interval, which is decided 

by the directional search. That is, if dH( p£) / d$\ 3 0, the next diagnostic interval is (p 2 , P2); 
if dH($) / d P2 a 0, the next diagnosis interval is (P 2 , P 2 ). 

8. Calculate the natural frequency with the mid-point stiffness p 2 of the new interval, in the 
procedure following from step 3 to step 7 until the following criterion for termination is satisfied 

IP^-P^I s e ( 58 ) 

where p^ denotes the nth diagnosed stiffness and pj"* 1 * the (n+l)th, from the next new interval 

following P^ . The value of e is chosen as 10 4 N in present example. Since the diagnosed P 2 

seems to be on the upper bound of the interval (see Table 2), we have to choose a larger 
reference value of p 2 , say P 12 . Let p^ 2 =2.5xl0 6 AT and 6 2 = 0-5xl0 6 AT. Therefore a second 

interval (2xl0 6 , 3xl0 6 )iV is created. 

9. Repeat the same procedure for the second model interval (2XK/ 5 , 3x1 O^N as for the interval 
( 0.5x10 6 , 0.8xl(f)N. The lower and upper bounds of diagnosed interval are also denoted by P 2 , 
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and 0*, respectively. 

10. Establish the convex model. Now one may notice that the diagnosed fc for the Sist interval 
tends to increase monotonically till O.SxlOW while for the second interval it decreased till 
2.0x10*//. This implies that the true (S ; may lie in the range between 0.8*1<?N and 2.0x10*//. 
It is straight-forward to determine the final possible small range of p 2 by employing the 
procedure indicated in steps 3 to 8. The above general procedure is illustrated in Fig. 2. 

Tables 2, 3, and 4 show the diagnosed intervals of every step for both initial hypothesized 
stiffness models and the final convex model, as well as the frequency distances at the 

corresponding mid-point values. 


Table 2 


step 

lower bound 
xlO- 6 

upper bound 
xlO- 6 

distance at 
midpoint 02 

distance at 
02 + A0 

sign of 
difference 

1 

0.5 

0.8 

1.0385 

1.0374 

' 

2 

0.65 

0.8 

0.8113 

0.8095 

' | 

3 

0.725 

— 

0.8 

0.6968 

0.6935 

1 ’ 1 

4 

0.7625 

0.8 

0.6403 

0.6391 

- 

5 

0.78125 

0.8 

0.6104 

0.6096 

- 

6 

0.79062 

0.8 
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Table 3 


step 

lower bound 
xKT 6 

upper bound 
xlO" 6 

distance at 
midpoint 

distance at 

Pi + ap 

sign of 
difference 


2.0 

2.5 

3.6642 


+ 

2 

2.0 

2.25 

3.1315 

3.1345 

+ 

3 


2.125 

2.8555 

32.8564 

+ 


2.0 

2.0625 

2.7120 

2.7154 

+ 

5 

2.0 

2.0313 

2.6414 

2.6514 

+ 

6 

2.0 

2.0156 

2.6135 

2.6158 

+ 

7 

2.0 

2.0078 

2.5915 

2.5929 

+ 


Table 4 


step 

lower bound 
xlO -6 

upper bound 
xlO -6 

distance at 
midpoint Pj 

distance at 
Pi+Ap 

sign of 
difference 


0.7906 

2.0078 

1.1006 

1.1043 

+ 

2 

0.7906 

1.3992 

0.2735 

0.2766 

+ 

3 

0.7906 

1.0949 

j 

0.1644 

0.1527 

- 

4 

0.9428 

1.0949 



+ 

5 

0.9428 

1.0188 

0.0518 

0.0491 

- 

6 


0.9998 



+ 
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From the above example, the following guidelines are suggested to construct the convex 

model: 

a) Initial Hypothesis : We first hypothesize one reference, nominal value of stiffness P 2 m 

and specify a certain 6 m . As a result, a interval as in Eq.(54), (P 2m * Pi,*"*” 

constructed and investigated by the present method. If the diagnosed p 2 belongs to this interval, 
the search process should identify it and thus the searching process should be terminated. This 
p 2 should be considered as the true value of stiffness. However, in most cases, due to the 

insufficient information, the reference value P 2 m may not be correct since the true value p 2 may 

not be included in the search interval. This situation will lead to the diagnosed value of stiffness 
p 2 to occur either at the upper or lower bound of the convex set. This would imply that more 
reasonable reference values have to be estimated and the previous convex model needs a 
correction or checking if the upper or lower bounds represent the true value. 

b) Correction: When the hypothesized convex model needs to be corrected, we have to 
create other convex model(s) which will include the true value of stiffness p 2 . Let us assume the 
diagnosed stiffness p 2 obtained in the initially hypothesized step lies at the lower bound of the 
interval, which means the true value of is equal or it is less than the lower bound of the 
interval. In the next step, we choose such a convex model among where all its values are less 
than those of the previous one. If the diagnosed stiffness p 2 of this convex model set also lies 
at the lower bound of the set, then we have to move the interval further to the left. Usually it 
is quite easy to create a convex model, such that the identified value of p 2 will suggest the upper 

bound of the interval, as long as the reference value P 2 „ is small enough. 
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c) Final Hypothesis: Once having such two alternative convex models, a final convex 
model can be constructed with the two previous diagnosed stiffnesses p 2 ’s representing the upper 
and lower bounds. As a result the "true" fi 2 is "caught" inside this model. 

2) Example 2 

In the second example, we will investigate the diagnosis of the alteration in the boundary 
condition of a square plate (a=b) with two adjacent edges clamped and the other two edges 
elastically supported (C-C-ES-ES) (see Fig.3). The uncertain rotational stiffnesses of the two 
elastically supported edges, 04 and |3 2 , are represented by the convex model, given in Eq.(l). As 
in the example 1 , the simulated situation (natural frequencies and corresponding underlying 
rotational stiffnesses) is constructed since experimental measurements have not yet been 
performed. 

Following the steps of imitating the real situation in the first example, we visualize that 
the set of measured nondimensional natural frequencies is given in Table 5. We note in passing 
that the frequency parameters in Table 5 were obtained by setting a 2 =lxl0 6 N, (5 2 =2xl0 6 N. The 
notations are same as that in example 1 . 
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Table 5 The first six nondimensional frequency parameters of the to-be-diagnosed plate, 
simulating the results of measurements 


Na_ 

28.678 

\ 2) 

62.882 

^3) 

94.693 

\«) 

117.162 

\s) 

148.056 

\s) 

191.336 


As discussed before, the nominal values of torsional stiffnesses are denoted by 
a 1} P x , a 2 , P 2 . Then the convex model set in Eq.(l) is as follows 


(59) 


C ,M = {a.,p,.,/ = l,2: 04=04=00, Pj = = |a 2 - =s b m , |P 2 - * bj 


or, simply 


C 2,„ = {a 2>p2 : l a 2 - a Z*l * 6 »> 102 - hm\ 


6 } 

m 


(60) 


since only 04 and p 2 , the stiffnesses of the elastically supported edges, are considered to be 
unknown. 

We will omit the intermediate steps to create the "final" convex model, since they are 
similar to those in example 1. Thus, we assume that we have arrived at the following convex 
model, o^O.lxlO 6 , 10xl0 6 )N, p 2 e(0.1xl0 6 , 10xl0 6 )N, which is large enough to include the 
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true values of stiffnesses a 2 and P 2 . The procedure, in detail, is as follows: 

1. We measure all the "necessary" natural frequencies of the actual plate (In this study, we 
assume that they are all provided and listed in Table 5.) 

2. We construct a convex set for altered boundary conditions. In this example it is 04, 
P 2 e(0.1xl0 6 , lOxlO 6 )//. The lower bounds of boundary stiffness intervals are denoted by a 2 and 

P 2 , whereas the upper bounds are indicated by and ($ 2 , respectively. 

3. We divide ( 0^, a 2 ) and (f} 2 , p 2 ) into a number of equal intervals (20 intervals in this study). 

For each interval’s separation point a 2 ( i = 1 , 2, 20) in (a 2 , a£) , we calculate the distances 

Hj between the measured and the calculated natural frequencies at all values of 
P 2 0 = 1, 2,—, 20), and select the combination (c 4, p£) at which the distance Hj reaches 

minimum. This step is actually equivalent to choosing (a£, (i£) among the 400 dividing points 

of ctj and (3 2 such that Hj attains its minimum. This may considerably save the storage due to 
the fact that some data becomes of no use once they are compared and need not to be stored any 
more. 

5. We redefine a searching rectangular range with the identified (c^*, p£) at ste P 4 as a center. 

The side length of this rectangular area is twice the divided interval length. It can be readily 
seen that the present area of a 2 and p 2 is one percent of the former one (see Fig. 4 & Fig. 5). 

6. We repeat the procedure of step 3 through step 5 until certain convergence criterion is 
satisfied. In this example, the criterion was chosen as the increment in values of 04 and p 2 
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between one search to another not to exceed lxl0 3 N. The numerical results for a’s and 
corresponding P’s at which the distances Hj reach minimum in different searching areas are listed 
in tables 6-9. 
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Table 6 The values of a and p at which Distances Hj reach Minimum (First Searching Region) 



a 

P at which H is min. 

H 

1 

0.100 

3.565 

2.3826 

2 

0.595 

2.080 

0.7270 

3 

1.090 

2.080 

0.4194 

4 

1.585 

1.585 

0.7438 

5 


1.090 

1.2431 

6 

2.575 

1.090 

1.7717 

7 

3.070 

0.595 

2.1220 

8 

3.565 

0.100 

2.3838 

9 

4.060 

0.100 

2.7720 

10 

4.555 

0.595 

3.2664 

11 

5.050 

0.100 

3.4329 

12 

5.545 

0.100 

3.7206 

13 

6.040 

0.100 

3.9796 

14 

6.535 

0.100 

4.2172 

15 

7.030 

0.100 

4.4345 

16 

7.525 

0.100 

4.6350 

17 

8.020 

0.100 

4.8179 

18 

8.515 

0.100 

4.9910 

19 

9.010 

0.100 

5.1491 

20 

9.505 

0.100 

5.2994 

21 

10.000 

0.100 

5.4382 
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Table 7 The values of a and p at which Distances Hj reach Minimum (Second Searching Region) 



a 

P at which H is min. 

H 

1 

0.5950 

2.4270 

0.5836 

2 

0.6445 

2.4270 

0.5006 

3 

0.6940 

2.3280 

0.4316 

4 

0.7435 

2.2290 

0.3654 

5 

0.7930 

2.2780 

0.2850 

6 

0.8425 

2.1790 

0.2144 

7 

0.8920 

2.1300 

0.1448 

8 

0.9415 

2.0310 

0.0834 

9 

0.9910 

2.0310 

0.0642 

10 

1.0410 

1.9810 

0.0569 

11 

0.1090 

1.9320 

0.1238 

12 

1 

1.1400 

1.8330 

0.1791 

13 

1.1890 

1.7830 

0.2424 



14 

1.2390 

1.7340 

liiiiifl 

15 

1.2880 

1.6840 

0.3651 

16 

1.3380 

1.6350 

0.4241 

17 

1.3870 

1.5850 

0.4824 

18 

1.4370 

1.5850 

0.5494 

19 

1.4860 

1.5850 

0.6155 

20 

1.5360 

1.5850 

0.6798 

21 

1.5850 

1.5850 

0.7438 
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Table 8 The values of a and |3 at which Distances Hj reach Minimum (Third Searching Region) 



a 

|3 at which H is min. 

H 

1 

0.9910 

2.011 

0.0119 

2 

0.9959 

2.006 

0.0068 

3 

1.001 

2.001 

0.0038 

4 

1.006 

1.996 

0.0083 

5 

1.011 

1.991 

0.0147 

6 

1.016 

1.981 

0.0212 

7 

1.021 

1.976 

0.0279 

8 

1.026 

1.971 

0.0341 

9 

1.031 

1.966 

0.0407 

10 

1.036 

1.956 

0.0466 

11 

1.041 

1.951 

0.0539 

12 

1.045 

1.946 

0.0598 

13 

1.050 

1.941 

0.0669 

14 

1.055 

1.932 

0.0723 

15 

1.060 

1.932 

0.0793 

16 

1.065 

1.932 

0.0873 

17 

1.070 

1.932 

0.0939 

18 

1.075 

1.932 

0.1009 

19 

1.080 

1.932 

0.1089 

20 

1.085 

1.932 

0.1159 

21 

1.090 

1.936 

0.1237 
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Table 9 The values of a and (3 at which Distances Hj reach Minimum (Fourth Searching Region) 



a 

P at which H is min. 

H 

1 

0.9959 

2.004 

■■ 

2 

0.9964 

2.004 

0.0041 

3 

0.9969 

2.003 

0.0032 

4 

0.9974 

2.003 

0.0029 

5 

0.9979 

2.003 

0.0022 

6 

0.9984 

2.001 

0.0021 

7 

0.9989 

2.001 

0.0006 

8 

0.9994 

2.001 

0.0008 

9 

0.9999 

2.000 

0.0005 

10 

1.000 

2.000 

0.0009 

11 

1.001 

2.000 

0.0016 

12 

1.001 

1.998 

0.0025 

13 

1.002 

1.999 

0.0027 

14 

1.002 

1.999 

0.0036 

15 

1.003 

1.997 

0.0038 

16 

1.003 

1.998 

0.0045 

17 

1.004 

1.996 

0.0051 

18 

1.004 

1.996 

0.0060 

19 

1.005 

1.996 

0.0064 

20 

1.005 

1.997 

0.0074 

21 

1.006 

1.996 

0.0083 
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6. DISCUSSION AND CONCLUSIONS 


The generalized Bolotin’s dynamic edge effect method, combined with the approach of 
the convex modeling of uncertainty, has been employed to estimate the convex models associated 
with the alteration of the boundary conditions of thin uniform rectangular plates. For 
determination of the natural frequencies the generalized version [7] of the Bolotin’s dynamic 
edge-effect method [6] was employed. The essential advantage of this general method is the 
possibility of finding the natural frequencies for any preselected pair of mode numbers. We have 
shown in the Appendix, by comparing the nondimensional frequency parameters with those due 
to Leissa [4] and Gorman [10], that the use of "postulate" (20) and of two Levy’s type solutions 
to approximate natural frequencies leads to very accurate results. 

The main idea implemented in the present paper is the convex modeling of uncertainty, 
as applied to the diagnosis of local modifications in the boundary conditions. In order to describe 
the alteration in the boundary stiffness of the plate, convex models are utilized to specify the 
possible realizations of boundary conditions. In the next step the multi-hypothesis decision is 
implemented, based on the minimum distance algorithm. 

For the sake of simplicity, plates with one or two edges elastically supported and the rest 
clamped were investigated. Two possible kinds of convex models for the stiffness of the 
elastically supported edge were constructed. The multi-hypothesis diagnosis, in conjunction with 
the generalized dynamic edge effect method, proved to be a fairly successful technique. 
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APPENDIX 


The numerical procedure of solving the first and the second auxiliary problems, resulting 
in Eq. (38) and Eq. (53), is carried out by applying Newton-Raphson method. Here Eqs. (36) 
and (51) are denoted as 


F(p, q) - 0 
G(p, q ) =0 


(Al) 

(A-2) 


They can be put in the expanded form as follows 


* 2 YiY 2 *i V 2 * 4 + 2n 2 y l y 2 s l S l p 2 R 2 * n\y 2 s 2 S 2 p 2 + ny J&sjpR* + 
iiYi * 2 C l s l pR 3 - ny^.c^pR 2 - ny^c^pR 2 * jiy^.C^pR + 
jiy l $ 1 C l s l pR ~ ny 2 $ 2 c x S 2 p - ny^^Sj + + $i+ 2 S i R + 

U 2 c?R - 2 UiC&R - U 2 S?R * U 2 C?R - = 0 


n 2 b l b 1 s 1 S 1 q 1 Z A + 2 ji 2 8 1 6 2 S 2 S 2 q 2 Z 2 + jfb&s&q 2 + nb^C^qZ 3 + 
nb^ 2 C 2 s 2 qZ 3 - JtS^c.S^Z 2 - nb^c^qZ 2 + nb^C&qZ + 
nb^ 2 C 2 s 2 qZ - nb^^q - n b,%c 2 S 2 q * ^%s 2 S 2 Z 2 + ^ t %s 2 Z * 
^,%c 2 Z - 2 y^ 2 c 2 C 2 Z - S}Z + t^C’Z - = 0 


According to Newton-Raphson method, the above two equations are approximately replaced by 


%?) - F. + (~)i (P ~ P) + (4~); (0 “ ^ = 0 
d/? 


(A-5) 


G(p, *) - G. + (^),. (p - P) + (fx (<7 “ <?i) - 0 


(A-6) 
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where p t and q t are the values of roots of Eqs. (A-l) and (A-2) in ith iteration and 
F. = F(p p <jr ; ), G. = G(p p q.). The (i+1 ) th iteration of the above solution yields 



n 

= n + 

a n G i 

~ a n^i 

(A- 7) 


i*+l 


a n a u 

~ a i2°21 




“ + 

a n a 22 

~ U ll G i 
- <3 n a 2i 

(A-8) 


where a u , a 12> a 21 , and a^ denote, respectively, 


— 

a = ( ) 

11 ' dp Jp “ Pl ’ qmq ‘ 

= ( dF ) 

a i 2 ^ dq pmp ‘ ’ q ' q ‘ 

(A-9) 

— 

- ( ^ G \ 

- ^~dp pm P* • qmq > 

^ fig • 9'9i 

(A- 10) 


The derivatives dF/dp, dF/dq, dG/dp, dG/dq were evaluated analytically using the computerized 
symbolic algebraic code REDUCE. In the p-q plane, the square defined by the lines p=m, 
p=m+l, and q=n, q=n+l, where m and n are any positive integers, there is only one root of Eqs. 
(A-l) and (A-2) (see Fig. 3). This allows us to find the natural frequency for any given pair of 
numbers m, n so long as the initial values of p and q are properly chosen. 

In Tables A-D the results are presented for the frequency parameters X=toa 2 qph/D 

computed using both the present generalization of Bolotin’s dynamic edge effect method, and the 
Rayleigh-Ritz method adopted by Leissa [4]. In cases where Leissa’s results are not available, 
the results by Gorman [10] are cited and denoted by the star. It should be noted that the mode 
sequence numbers, cited from Ref. [4], are somewhat different from those in Ref. [4], because 
repeated frequencies are considered here as one. 
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Table A 


Comparison of frequency parameters \=(x>a l 2 y ph/D (SS-C-SS-SS Plates) 


Mode 

Sequence 

a/b=2/3 

a/b=1.0 

a/b=1.5 

Ref. 3 

Present 



Ref. 3 

Present 

1 


15.578 

23.646 


42.528 

42.528 

2 


31.072 




69.003 

3 

44.564 

44.564 

58.646 

58.646 

116.267 

116.267 

4 


55.393 

86.135 

86.134 

120.996 

120.995 

5 

59.463 

59.463 

100.270 

100.270 

147.635 

147.635 

6 

83.438 

83.606 

113.228 

113.228 

184.101 

184.101 


Table B 


Comparison of frequency parameters X=coa J Jph/D (C-C-C-C Plates) 


Mode 

Sequence 

a/b=2/3 

a/b=1.0 

a/b=1.5 

Ref. 3 

Present 

Ref. 3 

Present 

Ref. 3 

Present 

1 



35.992 

35.113 

60.772 

59.809 

2 


41.191 

73.413 

72.899 

93.860 

92.682 

3 

66.143 

66.172 

108.27 

107.444 

148.82 

148.245 

4 

79.850 

78.766 

131.64 

131.701 

179.74 

179.978 

5 

100.85 





213.053 


i 


l After Gorman’s [10] result (*.=41.25) is multiplied by four, since the side lengths in Ref. 
10 are taken as 2a and 2b. 
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Table C 


Comparison of frequency parameters X=toa 2 \jph/D (C-C-C-SS Plates) 


Mode 

Sequence 

a/b=2/3 

a/b=1.0 

a/b=1.5 

Ref. 3 

Present 

Ref. 3 

Present 

Ref. 3 

Present 

1 

25.861 

25.692 

31.829 

31.438 

48.167 

47.622 

2 

38.102 

37.837 

63.347 

63.054 

85.507 

84.990 

3 

60.325 

60.111 

71.084 

70.879 

123.99 

123.657 

4 


65.447 

100.83 

100.460 

143.99 

143.681 

5 

77.563 

77.357 




211.054 

6 

92.154 

98.374 

130.37 

130.262 

214.78 

220.072 


Table D 


Comparison of Frequency Parameters X.=(off 2 <Jph/D (C-C-SS-SS plates) 


Mode 

Sequence 

a/b=2/3 

a/b=1.0 

a/b=1.5 

Ref. 3 

Present 

Ref. 3 

Present 

Ref. 3 

Present 

1 

19.952 

19.853 

27.056 

26.867 

44.893 

44.669 

2 

34.024 

33.907 

60.544 

60.549 

76.554 

76.290 

3 

54.370 

54.326 

92.865 

92.665 

122.33 

122.233 

4 


57.429 

114.57 

114.568 

129.41 

129.215 

5 

67.815 

67.694 

146.0* 

145.786 

152.58 

152.313 

6 

90.069 

90.013 

188.5* 

188.469 

202.66 

203.303 
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Clamped 



Fig. 1 A plate with three edges clamped and the fourth edge elastically supported (ES). 
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Fig. 2 The process of creating the "final" convex model. 
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CONVEX IDENTIFICATION OF BOUNDARY CONDITIONS 
BY FINITE ELEMENT METHOD 


Jianjie Fang and Isaac Elishakoff 


Center for Applied Stochastics Research 
and Department of Mechanical Engineering 
Florida Atlantic University 
Boca Raton, FL 33431-0991, USA 


Abstr act: The study deals with the identification of the boundary conditions by the finite element method. The 
identification of the convex model, to which the boundary stiffnesses belong, rather than the total reconstruction of the 
boundary conditions, is performed. Two example problems of the beams, one uniform and the other nonuniform, both 
Hampp/i at one end and elastically supported at the other, are considered and numerically evaluated in detail. 


1. INTRODUCTION 

The vibration problems of beams, either due to deterministic or random excitations, have 
been widely investigated. To the best of our knowledge, in most of these studies the boundary 
conditions at the ends of the beam are assumed as known (Gorman, 1975; Weaver, Timoshenko, 
and Young, 1990). Most investigators consider all possible combinations of boundary conditions 
at the ends, namely those of simply supported (SS), clamped (Q, free (F), or, more generally, 
elastically supported (ES) ends. However, in practical problems, engineers usually do not have 
the exact information associated with the boundary conditions; they may, however, possess some 
fragmentary knowledge about those boundary conditions, such as the probable ranges of the 
boundary springs’ stiffnesses. The natural question arises: How to identify the "true" boundary 
conditions based on the partial knowledge? 

Boundary conditions have traditionally been modeled as axial and torsional springs or 
translational and rotational springs at the boundaries. Identification of parameters such as spring 
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stiffnesses has usually been considered within the context of identifying mass and stiffness 
matrices from measured eigenvectors (Baruch, 1982; Kabe, 1985; Baruh and Meirovitch, 1985; 
Dong, Cheng, and Qi, 1991). Identification of eigenvalues and eigenvectors form the system 
response is a popular subject and there are several methods available, such as those in Refs. 
(Ibrahim, 1983; Ewins and Gleeson, 1982; Juang and Pappa, 1985). Baruh and Boka (1993) 
recently presented a procedure for the identification of the spring constants, by assuming that a 
number of eigenvalues and eigenvectors are known and making use of identified eigenvectors and 
orthogonality conditions. 

To the best of the authors’ knowledge, the convex treatment of the uncertainty in 
boundary conditions was addressed only by Ben-Haim and Natke (1992, 1993). They represented 
this uncertainty by convex models, which are sets of mathematical entities such as functions, 
vectors or matrices (for more details about convex modeling of uncertainty one may consult the 
monograph by Ben-Haim and Elishakoff (1990)). For example, in Ref. (Ben-Haim and Natke, 
1993), an adaptive diagnosis procedure was developed and furthermore, the performance of this 
procedure with uncertainty characterized by the convex models was studied. 

In this study, the incomplete or partial knowledge about the boundary conditions is 
characterized by the convex model. The main idea about convex models for uncertainty is stated 
in the following section. It must be stressed that diagnosis of the boundary conditions within the 
convex modeling does not imply the complete identification of the spring constants in the beam 
problem. Indeed, total reconstruction of the boundary stiffnesses is not only quite challenging 
but may be unnecessary for many operational or maintenance decisions. Rather, diagnosis of the 
boundary conditions may mean identification of the adversely affected boundary (x=0 or x-L, 
where x is the axial coordinate and L is the length of the beam), and moreover, estimation of 
limits on the magnitude of the stiffness change of the elastic strings modeling the actual boundary 
conditions. Once convex models of the boundary conditions are constructed, the diagnostic 
process in fact implies identifying, in the convex model, the element to which the actual 
boundary stiffness belongs. 

This will then allow us to precisely formulate the diagnosis as a discrete multi-hypothesis 
decision problem with attendant formulation of the adaptive termination of this algorithm. The 
boundary conditions are generally modeled in terms of both translational and rotational springs. 
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The present diagnosis is based on measurement of a certain number of natural frequencies of the 
beam. The integral part of the method is the subprocess of evaluating the frequencies of a beam 
with specified boundary conditions, which can be performed analytically for the uniform beams, 
or numerically for the non-uniform beams, by using the versatile finite element method or other 
approximate techniques. The class to which the "true" boundary stiffness belongs, can then be 
identified by the convex modeling and implementing the decision algorithm of the suitably 
chosen minimum distance between natural frequencies, or other suitable criteria. 


2. CONVEX MODELS OF UNCERTAINTY 
As indicated above, identification of the boundary condition means identifying in the 
convex model the element to which the actual boundary stiffness belongs. In this study, convex 
models developed by Ben-Haim and Elishakoff (1990) are employed to specify the degree of 
precision required in the diagnosis of the boundary translational and rotational stiffnesses. 

A convex model is a set possessing certain convexity properties. Each element of the set 
represents a possible realization of an uncertain quantity of interest. In our case, these elements 
are intervals within which the "true” stiffnesses, X^s, are possibly located, where i— 1,2 are 
associated with translational and rotational stiffnesses, respectively. Therefore, elements in this 
convex model can represent the acceptable uncertainty in the boundary condition. The diagnosis 
will be considered satisfactory when the uncertainty in the boundary condition has been reduced 
to an preselected level, i.e. when the interval with sufficient small length has been found to yield 
corresponding actual measured frequencies of the beam. This means that diagnosis of the 
boundary condition is in fact no more than identification of the interval to which the actual 
boundary stiffness belongs. It should be stressed again that the actual stiffness will not be 
identified; only the interval to which it belongs is sought by the diagnosis. 

While several different convex models are available (Ben-Haim and Elishakoff, 1990), the 

following simple model is particularly suitable for the beam problem. Let X/ * be the nominal 

stiffnesses, where m= 1, M is the total number of convex models, and i= 1, 2, corresponding 
to the translational and rotational stiffnesses, respectively. We define the convex model as a set 
of stiffness intervals whose centers are located at the nominal values. Specifically, the convex 
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model C is the following set of hypothesized stiffness intervals: 

C = (c m , m=l,2, —,AfJ 


( 1 ) 


where c m is 

c. -{*, : % -«■/■> K. i - 1,2} (2) 

In other words, c m is the set of translational and rotational stiffness intervals which deviate 
by as much as ±6 from the constant value K- m> . One recognizes that C is a uniform bound 
convex model. 


3. MULTI-HYPOTHESIS DECISION 

The reference values Kjf* in the convex models in Eq. (1) can assume any of the M 
different values Kf\ K l a) ...,K l m . Thus diagnosis of the boundary condition involves deciding 
which of the M convex models contains the true boundary stiffness. This decision is based on 
a multi-hypothesis formulation (Wald, 1947). 

In its simplest form the multi-hypothesis decision algorithm requires the choice of one 
representative or hypothesized stiffness from each element in the convex model. It should be 
noted that in its more general form, more than one hypothesized stiffness is chosen for each 
convex model. Let h m denote the hypothesized stiffness from the element c m . That is, h m E c m . 
Let o) w = (o> J w ,...w / w ) indicate a vector of the first J natural frequencies of a beam whose 
boundary stiffness is h m . Furthermore, let Q = (Q ; , ...Q;) represent the vector of measured 
natural frequencies. 

The "distance" from the measured natural frequency vector Q to the anticipated natural 
frequency vector to^ based on the mth hypothesis is : 

H m = &(Q, <o (m) ) ( 3 ) 

where denotes some measure of deviation between Q and u) m . One of the possible choices 
can be the following 
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^(Q,ioJ = |Q - oj ( ">I1 


(4) 


where ||jc |[ denotes a norm of vector x. 

It is reasonable to conclude that the true stiffness is likely to belong to the convex model 
whose anticipated natural frequencies are the "closest" to the measured ones. Let the index m 0 
satisfy: 


H - min H 
"• m " 


(5) 


The multi-hypothesis decision is that the true stiffness belongs to the element of c^. 


4. DIAGNOSIS PROCEDURE 

The decision algorithm implied by Eq. (5) can always reach a decision, no matter whether 
the least distance, H , is small or not. In other words, the multi-hypothesis decision will choose 

m a 

the most likely stiffness interval from the available options. 

As expected, in practical situation, it is assumed that all needed frequencies of the to-be- 
identified beam have been measured. We need to establish a convex model to which the 
boundary conditions belong. It should be pointed out that we have many definitions of the 
distance of frequencies given by Eq.(3), depending on the definition of the norm of vectors. The 
fact that the lower order natural frequencies can be measured more accurately than higher ones, 
allows us to reasonably define the distance as a sum of relative errors of frequencies over the 
first several natural frequencies, i.e, 


-E 



( 6 ) 


where N is the number of measured frequencies. Note that the distance in Eq. (6) does not 
represent the norm, yet it appears to be a preferable criterion, from the physical point of view. 
The diagnosis procedure is as follows: 
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1. Measure all the necessary natural frequencies of the beam. 

2. Construct the convex model of the unknown boundary conditions by taking advantage 
of the available preliminary knowledge of the range of values which can be reasonably expected 
to occur. In this study, a specific convex model with general form of Eq.(l) is established. We 
set 6j=0.005 N/m for the translational stiffness, and 6 2 =0.005 N-m for the rotational stiffness. 
This implies that the stiffness differing from each other by less than or equal to 0.005 N/m will 
be treated as coinciding with each other. However, the number of elements in the convex model, 
M, is taken as infinite, and as a result K/ m> and K 2 (mi can take infinite number of values. But they 
do have a finite range, which is .K i w E[0,10 4 ] EIJL 3 N/m, -flT/^EfO.lOO] EIJL N-m in our example 
problem, where [a,b] denotes an interval. Since there are infinite numbers of elements, we can 
not afford testing all the elements in the convex model. The adaptive process of producing sub- 
intervals from original large intervals is utilized to reach the element of interest efficiently, 
without the necessity of checking all the elements. In this process the lower bounds of the 
boundary stiffness intervals are denoted by K.* and K 2 a , and the upper ones by K t and K 2 
respectively. We note that in the first procedure K“= 0 N/m, K b - 10 4 EIJL 3 N/m, K 2 = 0 N-m, 
^=100 EIJL N-m. 

3. Divide intervals [KJ, K ,*] and [K 2 °, K^\ into a number of equal intervals (20 intervals 
were chosen in this study). 

4. Compute the distances & between the measured frequencies and the ones 
corresponding to hypothesized stiffnesses according to the general definition given in Eq. (3) for 
every interval separation point KJ (i- 1,2,..., 20) in [j K", Kf\ at all Kj’s (j= 1,2,. ..,20) in \K 2 , K 2 ], 
by employing the finite element method (Petyt, 1990). 

5. Identify the point (K{\ K 2 ) from the multi-hypothesis decision of Eq. (5) after 
obtaining all the distances in step 4. 

6. Reset a searching range of K 2 and K 2 , which is a rectangle in two dimensional space, 

with the identified vector (<’, k£) from step 4 as the center, and the side length of this 

rectangular is taken as twice the divided interval length of the preceding step. It can be readily 
seen that the present area of K ; and K 2 is one percent of the preceding one. 
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7. Repeat steps 3 through 6 until an interval is found such that, on the basis of minimum 
distance algorithm decision Eq. (5), H m is smaller than a specified threshold value 0. In other 
words, the process of diagnosis is terminated if the following inequality is satisfied 

H C 7 ) 

m 

In the present example, this step is implemented as choosing (K^, K 2 ) among the 400 dividing 

points of Kj and K 2 such that H is minimum. 

As we can see, one of the cornerstones of the method is the ability to determine the 
natural frequencies of beams under arbitrary boundary conditions. The finite element method, 
which is addressed in Appendix A, is applied to solve this problem. 


5. NUMERICAL EXAMPLE 

As an example, we will investigate the beam which is clamped at one end and elastically 
supported at the other. 

5.1 Uniform Beam 

Consider a uniform beam with length L, Young’s modulus E, material density per unit 
length p, moment of inertia I and cross-sectional area A. In this case, we visualize that only the 
first two nondimensional natural frequencies were measured, which are given as follows 

\ = 18.0 (8) 
\ = 50.0 


The relation between the natural frequencies and their nondimensional counterparts A. ; is 


0l> = 

i 


\ EI_ 
TN pA 


(9) 


where to, is the natural frequency of the beam. 

We note in passing that the exact frequencies for the uniform C-C (i.e., with both ends 

clamped) beam are 
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( 10 ) 


\ = 22.37 \ = 61.67 

\ = 120.9 \ - (i+ 1/2) V , /or i*4 


and for a C-F uniform beam are 

X t = 3.516 \ = 22.03 

Xj = 61.70 X. - (i - 1/2) 2 * 2 , /or is 4 


( 11 ) 


These cases are of interest since C-C and C-F are two extreme cases of C-ES. 

The results of the identified intervals in various steps are summarized in Table 1. 

Table 1 


step number 

Interval for k 3 

Interval k 4 

1 

[0, 10]xl0 3 

[0, 100] 

2 

[1.0,2-OJxlO 3 

[1,10] 

3 

[I.l5,1.25]xl0 3 

[4.5, 5.5] 

4 

[1.22,1.23]xl0 3 

[4.9, 5.0] 

5 

[1.2215, 1.2225]xl0 3 

[4.965,4.975] 

6 

[1.2216, 1.2217]xl0 3 

[4.970,4.972] 


where the nondimensional stiffnesses are 


1 EL 


K = 


K,L 

El' 


k, = 


K 3 L 3 

~el' 


k, = 


K<L 

EL 


( 12 ) 


Thus, the final result is k 3 E.[ 1.2216, 1.2217]xl0 3 , ^E[4.970, 4.972]. One can take the 
mid-points of these intervals k 3 =1.22165xl0 3 , k 4 =4.971 as the identified boundary stiffnesses. 
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The first three nondimensional natural frequencies, obtained by finite element method for above 
k 3 and k 4 , are ^=18.0000, ^=49.9999, X 3 =94.0661. 

The above results are verified by the weighted residuals method which is given at 
Appendix B of the paper. The single term weighted residuals method yields X 2 =18.0637, which 
constitutes an error of only 0.3539%. This is also verifiable through the exact solution, which 
is obtained by solving the following characteristic equation (Weaver, Timoshenko, and Young, 
1990): 

0 1 0 1 

1 0 1 o 

-X T cosy/X - & 3 siny/X X T siny/X-k,cosy/X X 7 coshy/X - kjSinhy/X X 7 sinhy/X - fc 3 coslVx 

-VXsim/X + Jk 4 cosy/X -V^cosy/X - A: 4 sini/X Jki sinhy/X + fc 4 coshy/X y/Xcoshy/X + jfc 4 sinhy/X 

= -ifc 3 Jt 4 cos 2 y/X - X 2 cos 2 y/X + 2fc 3 £ 4 cos y/x cosh y/x - 2X 2 cos y/X cosh y/x - k 3 k 4 cosh 2 y/X 
- X 2 cosh 2 y/X - 2 k 3 y/X cosh y/X sin y/X - 2X 4 X y/X cosh y/X sin y/X - ^ 3 X 4 sin 2 i/X - X 2 sm 2 ,/X 
+ 2X 3 y/X cos ^X sinh y/X - 2k A X y/X cos y/X sinh y/xT + fc 3 k 4 sinh 2 y/X + X 2 sinh 2 y/X = 0 ( 13 ) 

Substitution of £,=1.22165xl0 3 , fc 4 =4.971 into Eq. (13) gives the following first three non- 
dimensional frequencies: X i= 17.9996, X 2 =49.9913, X 3 =94.0085, which are extremely close to the 
measured ones. 

5.2 Non-Uniform Beam 

Consider a non-uniform beam with length L, Young’s modulus E, the moment of inertia 
I and cross-sectional area A are varying proportionally along the axis ox, with I, and A, at one 
end and I 2 and A 2 at the other. Similar to the case of uniform beam, we visualize that only the 
first two nondimensional natural frequencies are measured, which again are given in Eq. (8). For 
the special case A 2 =A,/2=A/2, and I 2 =I,/2=I/2, the similar nondimensional stifftiess intervals are 
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listed as follows 


Table 2 


step number 

Interval for k 3 

Interval k< 

1 

[0, 10]xl0 3 

[0, 100] 

2 

[0,1.0]xl0 3 

[1,10] 

3 

[0.55,0.65]xl0 3 

P.3] 

4 

[0.640, 0.650JX10 3 

[2.35,2.45] 

5 

[0.6400, 0.6410]xl0 3 

[2.410,2.420] 


Thus, the final result is £ 3 E[0.6400, 0.6410]xl0 3 , £*£[2.410, 2.420]. The first three 
nondimensional natural frequencies corresponding to mid-intervals, i. e. to £ 3 =0.6405xl0 3 , 
£.,=2.415, are ^=17.9993, ^=50.0013, X, =94.3288. 


6. DISCUSSION 

As one can expect, the availability of the information, namely of the number of the 
measured natural frequencies, will certainly have an effect on the identification process of the 
boundary conditions. It appears to be instructive to investigate the effect of adding the 
information of higher order frequencies to the criterion of frequency distance. The following two 
cases will be investigated. If the information on the third nondimensional frequency X i= 94.0 is 
added to the criterion of frequency distance Eq.(6), i. e. for case of measured natural frequencies 
being equal to X 7 =18.0, ^=50.0, X 3 =94.0, the results of the identified intervals in various steps 

are listed in Table 3. 
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Table 3 


step number 

Interval for k 3 

Interval k 4 

1 

[0, 10]xi0 3 

[0, 100] 

2 

[1. 0,2.0] xlO 3 

[1.10] 

3 

[1.15,1.25]xl0 3 

[4.5, 5.5] 

4 

[1.21,1.22]xl0 3 

[4.95,5.05] 

5 

[1.215, 1.216]xl0 3 

[4.975,4.985] 

6 

[1.2153, 1.2154JX10 3 

[4.977,4.979] 


where k 3 and k 4 are the nondimensional stiffnesses given by Eq. (12). For the mid-interval 
stiffnesses £ 3 =1.21535xl0 3 and £<=4.978, the first three theoretical nondimensional frequencies 
are ^=18.0000, X2=49.9834, ^=94.0001 respectively. 

Let the third nondimensional frequency be taken as X 3 =94.0661, which is the actual 
nondimensional natural frequencies corresponding to the "identified" boundary stiffnesses 
£ 3 =1.22165xl0 3 , £<=4.971. The results of the identified intervals in various steps are as follows: 

Table 4 


step number 

Interval for k 3 

Interval k 4 

1 

[0, 10]xl0 3 

[0, 100] 

2 

[1.0,2.0]xl0 3 

[1,10] 

3 

[I.l5,1.25]xl0 3 

[4.5, 5.5] 

4 

[1.215, 1.225]xl0 3 

[4.95,5.05] 

5 

[1.2210,1. 2220]xl0 3 

[4.965,4.975] 

6 

[1.2216, 1.2217]xl0 3 

[4.970,4.971] 
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where k 3 and k 4 are the nondimensional stiffnesses given by Eq. (12). Thus the process 
converges to the same final interval as it should be. 

Assume now that the measurement shows the same first two natural frequencies but with 
different third frequency, namely X 3 =95.0. This simulates the imprecision in the measurements. 
The results of the convex identification leads to intervals in various steps as follows: 


Table 5 


step number 

Interval for k 3 

Interval k 4 

1 

[0, 10]xl0 3 

[0, 100] 

2 

[1.0,2.0]xl0 3 

[1.10] 

3 

[1.25,1.35]xl0 3 

[4.5, 5.5] 

4 

[1.31,1.32]xl0 3 

[4.8, 4.9] 

5 

[1.316,1.317]xl0 3 

[4.87,4.88] 

6 

[1.3167, 1.3168]xl0 3 

[4.875,4.876] 


where k 3 and k 4 are the nondimensional stiffnesses given by Eq. (12). If one decides to let the 
"identified" frequency to be at the middle of the intervals, i.e. k 3 = 1-3 1675x10 and k 4 =4.8755, the 
first three corresponding theoretical frequencies will be X ; = 18.0000, X. 2 =50.2287, X 3 =95.0002. 
As is seen, the values of k 3 has changed by 8%, and values of k 4 has been changed by 2%, 
although the error in the third frequency constituted slightly above 1%. This shows that for 
accurate identification of the boundary conditions one needs much precision in the measured 
natural frequencies. 

The above results show that the identified stiffnesses will change with the additional 
information on the measured higher order frequencies added to the criterion of the frequency 
distance Eq. (6). 

To understand the present convex identification method, one may get some additional 
insight by contrasting it with the deterministic identification. By the finite element method, for 
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the known stiffnesses k 3 and k 4 , one can always numerically obtain the corresponding 
nondimensional frequencies. Fig. 1 shows the variation of the first nondimensional frequency 
\ for different values of k 3 and k 4 . Similarly, we can obtain frequencies X* X 3 , -, X, for any 
n for different values of k 3 and k 4 . By making cuts at X ; =18 and X 2 =50 in the two three- 
dimensional figures for X 2 and X 2 , we obtain two curves: curve (a) gives the combination of 
values k 3 and k 4 for which the first frequency X 2 is constant, namely equals 18. Curve (b) is 
associated with the second frequency T^fk# k 4 )=50. The intersection point of the two curves in 
Fig. 2 must identify the values of k 3 and k 4 if there are no errors in measurements. However, it 
is important to keep in mind that these two curves in Fig. 2 are not extremely accurate because 
in figures of X 2 and X 2 the nondimensional frequencies X 2 and X^ are not continuously changing 
with k 3 and k 4 , and some error is present. Furthermore, to reach high accuracy, extremely large 
amount of points are needed to obtain a collection of k 3 and k 4 for which X 2 =const. and X 2 =const., 
i.e. equal frequency curves. On the other hand, when higher order natural frequencies are 
measured, the deterministic identification method may be unapplicable. The reason is that the 
equal frequency curves obtained through higher order frequencies may not intersect at the same 
point where the curves X 2 =const. and X 2 =const. intersects due to the unavoidable measurement 
errors. However, the convex identification is still able to treat this effect. This is the advantage 
of convex modeling since it takes into account the uncertainty in the measurements. 

7. CONCLUSION 

In this problem, a special kind of convex model of boundary conditions of a beam is 
utilized to specify their possible realizations. The diagnostic process in fact implies identification 
of the interval to which the "true" boundary conditions belong. This process is realized by 
implementing the multi-hypothesis decision technique which is based on the minimum distance 
algorithm. The beams with one end clamped and the other end elastically supported are 
investigated. We assume that a certain number of natural frequencies of the beam were measured 
experimentally. The evaluation of the frequencies of a beam with specified boundary conditions 
is performed numerically by the finite element method. The final identified boundary stiffness 
intervals are checked by comparing the calculated natural frequencies with the measured ones. 


2-13 


It is found that they are extremely close to each other. The effect of the additional information, 
namely the number of the measured natural frequencies, on the identification results, is also 
investigated in this study. The present convex model identification based on multi-hypothesis 
decision can be utilized in identifying damaged regions in the boundary conditions in plates and 
shells. In addition to versatile finite element techniques, one can use some analytical techniques 
for uniform structures. Such a method of calculation of natural frequencies in rectangular plates 
with step-wise discontinuities in rotational stiffness along the edges was developed by Gorman 
(1993) at the suggestion of one of the present writers. 


ACKNOWLEDGMENT 

This study has been supported by the NASA Kennedy Space Center, through Cooperative 
Agreement No. NCC10-0005, SI to the Florida Atlantic University; technical monitor: Mr. R. 
Caimi. Authors appreciate fruitful discussion with Professor Y. Ben-Haim of the Technion- Israel 
Institute of Technology. 


2-14 


REFERENCES 


Baruch, M., 1982, "Correction of Stiffness Matrix Using Vibration Tests," AIAA Journal, Vol. 
20, No. 3, pp441-442. 

Baruh, H. and Meirovitch, L., 1985, "Parameter Identification on Distributed Systems," Journal 
of Sound and Vibration , Vol. 101, No. 4, pp551-564. 

Baruh, H. and Boka, J.B., 1993, "Modeling and Identification of Boundary Conditions in Flexible 
Structures," International Journal of Analytical and Experimental Modal Analysis, Vol. 8, 
ppl07-117. 

Ben-Haim, Y. and Elishakoff, I., 1990, Convex Models of Uncertainty in Applied Mechanics, 
Elsevier Science Publishers, Amsterdam. 

Ben-Haim, Y. and Natke, H.G., 1992, "Diagnosis of Changes in Elastic Boundary Conditions in 
Beams by Adaptive Vibration Testing," Archives of Applied Mechanics, Vol. 62, pp210-221. 
Ben-Haim, Y. and Natke, H.G., 1993, "Sequential Adaptation in Estimating Elastic Boundary- 
Condition Influence Matrices," Journal of Dynamic Systems, Measurements and Control, 

Vol. 115, No.3, pp370-378. 

Cheung, Y.K., 1993, Finite Element Methods in Dynamics, Kluwer Academic Publications, 
Dordrecht. 

Dong, X.Q., Cheng, Y. and Qi, J., 1991, "Identification of Structural Boundary Condition," 
Proceeding of 9th International Conference of Modal Analysis, Firenze (Florence), Italy, pp993- 
998. 

Elishakoff, I. and Fang, J., "Diagnosis of Local Modifications in the Boundary Conditions of a 
Rectangular Plate," Journal of Computer Methods in Applied Mechanics and Engineering, 
(submitted). 

Ewins, D.J. and Gleeson, P.T., 1982, "A Method for Modal Identification in Lightly Damped 
Structures," Journal of Sound and Vibration, Vol. 84, No. 1, pp57-59. 

Gorman, D.J., 1975, Free Vibration Analysis of Beams and Shafts, Wiley, New York. 

Gorman, D.J., 1993, "Free Vibration Analysis of Rectangular Plates with Step-Wise 
Discontinuities in Rotational Stiffness Along the Edges," Journal of Sound and Vibration, Vol. 


2-15 


166, pp397-408. 

Ibrahim, S.R., 1983, "Computation of Normal Modes from Identified Complex Modes," AIAA 
Journal, Vol. 21, No. 3, pp446-451. 

Juang, J.N. and Pappa, R.S., 1985, "An Eigensystem Realization Algorithm for Modal Parameter 
Identification and Model Reduction," Journal of Guidance Control Dynamics, Vol. 8, NO. 5, 
pp62G-627. 

Kabe, A.M., 1985, "Stiffness Matrix Adjustment Using Mode Data," AIAA Journal, Vol. 23, No. 
9, ppl431-1436. 

Leissa, 1969, A.W., "Vibrations of Plates," NASA SP 160, 1969. 

Petyt, M., 1990, Introduction to Finite Element Vibration Analysis, Cambridge University Press. 
Wald, A., 1947, Sequential Analysis, John Wiley. Re-issued by Dover Press, 1973. 

Weaver, W., 1986, Structural Dynamics by Finite Elements, Englewood Cliffs, Prentice-Hall. 
Weaver, W. Jr., Timoshenko, S.T. and Young, D.H., 1990, Vibration Problems in Engineering, 
John Wiley and Sons, New York. 


2-16 


APPENDIX A. VIBRATION OF ELASTICALLY SUPPORTED BEAMS 
BY FINITE ELEMENT METHOD 

In order to perform the convex identification procedure, one needs to have an analytical 
or numerical algorithm for evaluation of natural frequencies of the structure under consideration. 
In our previous study (Elishakoff and Fang), devoted to elastically supported plates, the 
generalized Bolotin’s dynamic edge effect method was utilized. For plates, one can use one of 
the numerous methods reviewed in the monograph by Leissa (1969), or employ the analytical 
approach based on the principle of superposition, by Gorman (1993). In this section the versatile 
finite element method will be discussed for vibration frequency evaluation. Whereas the finite 
element method for beams is exposed in several books, (Cheung, 1993; Petyt,1990; Weaver, 
1986; Weaver, Timoshenko, and Young, 1990), none of them unfortunately deal with elastically 
supported beams. This gap is closed in this paper. 

Strain energy of elastic deformations stored in an element of length 2a reads 




(A.1) 


where E is the modulus of elasticity, and I t is the second moment of area of the cross-section 
about the z-axis. The kinetic energy of the element is 


T - 1 l‘pA* 2 dx 


(A.2) 


where p is the mass density, A is the cross-sectional area, dot denotes the derivative with respect 
to time. The element shown in the figure has a total of four degrees of freedom, displacements 
and rotations at each end of the beam (Petyt, 1990). The displacement function can thus be 
represented by a polynomial having four constants, namely 

(A.3) 


v = + 0 ^% +a 3 | 2 + a 4 | 3 , | = x/a 

This expression can be written in the following matrix form 
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V ' [Pd)] r (al 


(A.4) 


where 

Wl)] r - [i ? V I 3 ] 

[af = [a, a, a, aj 

Differentiating (A.3) yields 

a0 = a— - — = ctj +203^ + 3 a 4 ^ 2 
2 dx d% ^ 

Evaluating (A.3) and (A.6) at ^=-1, and ^=1 results in 


’ V 1 ' 


1 

-1 

1 

-1 


a i 

00 


0 

1 

-2 

3 


a z 

V 2 


1 

1 

1 

1 


«3 

00 


0 

1 

2 

3 


a 4 

Z 2_ 








where v 1= v(-l), 0 zl =0(-l), v 2 =v(l), ©, 2 = 0 ( 1 ), or, in matrix form 

{v) = [A] e (a) 

where 

{v) = (v i; ad Zi , v 2 , ad^F 

Solving (a) from Eq. (A.8) leads to 

(a) = [A], 1 (v) 

where 


(A. 5) 


(A.6) 


(A.7) 


(A.8) 


(A.9) 
(A. 10) 
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2 12-1 
\ -3 -1 3 -1 

4 0-10 1 

11-11 


Eq. (A. 10) can be written in the different form 


{«) = [c] e H 


where 


{vf = 


V, 0 V. 0 

1 z, 2 z. 


and 


M. - 4 


2 a 2 -a\ 

1-3 -a 3 -a | 

0 -a 0 a 

1 a -1 a 


Substituting (A12) into (A. 4) yields 


v - wnmH 

Eq. (A. 15) can be expressed in the form 

V = [«(!=) ]{v>. 


where 


[N($] = [N&) aN 2 (%) N 3 m) aNJ§] 


The displacement functions in (A. 17) are given by 


(A. 11) 

(A. 12) 
(A. 13) 

(A. 14) 

(A 15) 
(A. 16) 
(A. 17) 
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*,(l) -I(2-3$-£) 
NJ&) - 1(1 -I -I 2 *? 5 ) 
IV, (0 -i(2 *3? -|>) 
JV,® -i(-l -1*0*0) 


(A. 18) 


Substituting the displacement expression (A. 16) into the kinetic energy (A.2) results in 

T e = i J B pAv 2 (ix = .1. J"°pA v 2 a<#j 
= I<v£pa fA(|) [*(1) ] r [tf(|)]di{v}, (A. 19) 

2 J-i 

= j <*? [m] e W 


Therefore, the element’s inertia matrix is given by 

k - p « / ) m i m r t *(?) i ^ 


(A. 20) 


Substitution of the displacement expression (A. 16) into the strain energy (A.1) results in 


U = - [ °EI {—fdx = 1 [ l EI z ±.^fa 

• 2 J - 2 dx 2 2 J-i a 4 d £ 2 


ii v 

dx 


1 r 1 1 , dV 






= I Mi JL j ; /j® t at- " © r [ M, 


(A.21) 


- |<v£w>>. 


The element stiffness matrix is 
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(A.22) 


[*L - 4 [""(?) 

a i J * l 


For a uniform beam, the following element inertia matrix can be obtained from Eq.(20) 





78 

22a 

27 

-13a 

pAa 

22a 

8a 2 

13a 

-6a 2 

105 

27 

13a 

78 

-22a 


-13a 

-6a 2 

-22a 

8a 2 


(A.23) 


The element stiffness matrix 


is obtained from Eq. (22) 


[*], - 



3 

3a 

-3 

3a 

El 

3 a 

4a 2 

-3a 

2a 2 

2 a 3 

-3 

-3a 

3 

-3a 


3a 

2a 2 

-3a 

4a 2 


(A. 24) 


For a non-uniform beam with linearly changing cross-sectional area and inertia of moments, the 
element inertia matrix is obtained through the use of the computerized algebraic code 
MATHEMATICA 


[»L 



‘ 78 

22a 

27 

-13a 


-42 

-8a 

0 

a 

^1 +J ^2 P <3 

22a 

00 

fc 

Ni 

13a 

-6a 2 

^ 44 2 -A 1 p a 

-8a 

-2a 2 

a 

(A- 

2 105 

27 

13a 

78 

-22a 

2 105 

0 

a 

42 

-8a 


-13a 

-6a 2 

-22a 

8a 2 


a 

0 

-8a 

i 

rs 

<N 


The element stiffness matrix is 


[U 



3 

3a 

-3 

3a 


0 

-a 

0 

a 

A +/ 2 £ 

3a 

4a 2 

-3a 

2a 2 

^ A “A £ 

-a 

-2a 

a 

0 

2 2a 3 

-3 

-3a 

3 

-3a 

2 2a 3 

0 

a 

0 

-a 


3a 

2a 2 

-3a 

4a 2 


a 

0 

-a 

2a 2 


(A.26) 
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where A lf A 2 and / ; , I 2 are cross-sectional areas and moments of inertia at the two ends of the 
element, respectively. 

Let us consider the expression of the strain energy of the elastic supports. The whole 
beam is shown in Fig. 3. For element 1, the strain energy of the elastic support is 


V? - A *.* 1 


(A.27) 


But 

v l = [1 0 0 0]M e 
0 2 = [ 0 1 0 0] M e 

Therefore 



i 


10 0 0 

< 

- M 

II 

0 

0 

[ i o o o ] M, = MI 

0 0 0 0 
0 0 0 0 


0 


0 0 0 0 


0 


0 0 0 0 

e z 2 , = WI 

1 

0 

[o l o o ] M c = Ml 

0 10 0 
0 0 0 0 


0 


0 0 0 0 


{v) 


M_ 


(A.28) 


(A.29) 


(A.30) 


Therefore the stiffness matrix of elastic support at end 1 is 



^000 
0 k 2 0 0 
0 0 0 0 
0 0 0 0 


(A. 31) 
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Similarly the stiffness matrix of elastic support at end 2 is 




(* + i) _ 


0 0 0 0 

0 0 0 0 

0 0 k 3 0 

0 0 0 k A 


(A. 32) 


Note that there are no support contribution at the elements 2 through N-l. Therefore the stiffness 
matrix at element 1 is 


id) 


and that at element N is 


[*£ - [*L + [*]i 


[kfe = + 


(A.33) 


(A.34) 


i =2, 3, -,N - 1 


(A.35) 


And for element 2 through N-l, 

[*]i - [*L> 

Following the general procedure of assembling the elemental mass and stiffness matrices 
into global mass and stiffness matrices, we can obtain [M] and [X]. The eigenvalue problem is 

[tf-co 2 M]{3>} = 0 ( A36 ) 

where to is the natural frequency of the beam. The nondimensional natural frequencies can be 
obtained by 


X. = to L 


pA 

El 


(A. 37) 


In the present study, subspace iteration method is applied to solve eigenvalue problem. 
For the case of two element approximation, the mass matrix reads 
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‘ 78 

22a 

27 

-13a 

0 

0 



22a 

00 

13a 

-6a 2 

0 

0 

— 

r m i _ 

W - ‘ 105 

27 

-13a 

13a 

-6a 2 

156 

0 

0 

16a 2 

27 

13a 

-13a 

-6a 2 




0 

0 

27 

13a 

78 

-22a 
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0 

-13a 

-6a 2 

-22a 

8a 2 


The element stiffness matrix is 
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3 +k{ 

3a 

-3 

3a 

0 

0 



3a 

(4 + k^)a 2 

-3a 

2a 2 

0 

0 

— 

r , , El 

-3 

-3a 

6 

0 

-3 

3a 


[k] = —— 
L J< 2a 3 

3a 

2a 2 

0 

8a 2 

-3a 

2a 2 



0 

0 

-3 

-3a 

3 +*/ 

-3a 



0 

0 

3a 

2a 2 

-3a 

(4+fc 4 ')a 


where k', k[, ki, k' are as follows 


— 

k( = k x 

1 El 

k[ = k — 
2 El 

— 

k' = k 3 

3 El 

ii 

Sis 


(A. 38) 


(A.39) 


(A.40) 
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APPENDIX B. WEIGHTED RESIDUALS METHOD 
The governing differential equation of a beam reads 


d 2 

dx 2 


t 


\ 


Em 


d 2 ™ 

dx‘ 


=« 

dt 2 


(B.l) 


where E is Young’s modulus, I(x) is the moment of inertia, p is the material density, A(x) is the 
cross-sectional area, w(x,t) is the transverse displacement, x is the axial coordinate and t is time. 
The boundary conditions are 


„ u 




= 0 


„ dw -T! \ d 2 ™ n 


„ d ( pt/ \ d 2 ™ 

K,w-i 3 —Em— 


\ 


= 0 




at x =0 


at x =0 


at x - l 


at x = l 


(B.2) 


Here, K, and K } are the stiffnesses of the translational springs, and K 2 and K 4 are the stiffnesses 
of the rotational springs. In addition, artificial parameters, y 1} y* Y* Y* taking value zero or 
unity, are introduced to model all possible idealized boundary conditions. 

For a nonuniform beam with linearly changing cross-sectional area and inertia of 

moments, we have 


A(x) =Aj 

m = a 






(B.3) 


where 
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a = 


A 2 _ E 

Ai Iy 


(B.4) 


Note that a=l corresponds to the case of a uniform beam. 

We introduce the nondimensional space coordinate | and time coordinate x: 




PAi 


(B-5) 


Therefore the differential equation (B.l) becomes 




(B.6) 


Then the boundary conditions can be rewritten as 


-° " t-0 

ai a | 2 

*,i^- Yl [l-(a-l)|]^ =0 at |=0 

t ,».- Y ,±{[ l *( a - l )|]^} -0 «»?- 1 

ai ai 2 


* 4 ±l n< [l.(a-l)|]i^ =0 
4 ai 4 ai 2 


at i = 1 


(B.7) 


where the nondimensional stiffnesses are 


1 EL 


k ~^ L K - K ' L ' 

2 El* * El, 


k. = 


K<L 

EL 


(B.8) 


To obtain some insight into the problem, a single-term weighted residuals approximation 
associated with the computerized symbolic algebra is used in this study. We first construct a 
mode-shape function which satisfies both geometric and dynamic boundary conditions in the form 
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w(i) = c 0+ c 1 i + c 2 i 2 + c J i 3 + i 4 


(B.9) 


Coefficients c,’s in the mode shape function, obtained through the use of the computerized 
symbolic algebraic code MATHEMATICA, are given as follows 





(B.10) 


where 


A„ = -(-24 Yl Y3* 2 +48a Yl Y 3 A: 2 -2Aa\y 3 k 2 - Uyyjc, -Syjcjc, + 2ay l k 7 k 3 )(12ay 4 +4 kj 
- (24a YlY4 A: 2 - 12 a\y& + 12 YlY ^ 4 + 18 Yl ^ 4 - 6a Yl ^ 4 )(12 Y3 -36a Y 3 + kj 
A l = -(-Uyj^ +240^3^ + 12 YlY2 ^ 3 -2y 7 k l k)(l2ay 4 +4 k} 

-< 12 Y 3 - 36a Y 3 + ^ 3 )(12a Y2 Y 4 A: 1 + 6y^k^) 

A 2 = -<-6y 3 k t k 2 + 12a Y3 ^ 2 + 6y l k 2 k 3 - k l k 2 k 3 )(12ay 4 + 4^J -(12y 3 - 36a Y3 + k^ayjc^ + 3*^ 
A 3 = -{2y 3 kjc 2 - 2ay 3 k l k 2 + 2^ + 2 Yl ^ 3 - 2a Yl ^ 3 + k 3 k^{12cty 4 + 4^ 

-(12 Y 3-36a Y3 +k 3 )(-2ay 4 k l k 2 -2y 7 k l k 4 -2kjcjc^ 

A = \2a\y 4 k 3 k 2 + 12a Y2 Y 4 ^ 3 + 24a Yl Y 4 ^ 3 - 12a +4a Y4 ^3 - 12 y 2 y 3 ^ 4 + 240^^ 

-6Y 3 ^ 4 +18a Y 3^ 4 + 12 YlY ^ 4 +*(£££< + 18 Y M*4 -6ayjc^k 4 + kjcjc 3 k 4 

(B.ll) 

For the case of uniform beam the coefficients A; read 
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A 0 = 12 Yl {[(5^ 3 + 24 Yj)y 4 *(k 3 +24 y3 )k 4 ]k 2 + 3y 2 [A + 8 Y3 )* 4 + 4 Y / 3 ]} 

A t = 2 y 2 {[(^ +48y^)k 4 + 6(* 3 + 12^]^ - 24{* 4 - 3^^) 

K 

4 - 2 4 

2 T 1 (B.12) 

A 3 = -2{[(5^3 + 24 Y3 )y 4 + (k 3 + 24 Y3 )* 4 ]*A + 3[(* 3 + 8 Y3 )* 4 + 4 Y A] Y A> 

A = {[(*, + 12 Y j)^ 4 + 4 (*, + 3 Ys ) Y JA: 1 + 12(* 4 + 

+4{[A +3^ +3 Y aA + 3 Y A^ 4 hr 2 


We seek a solution of Eq. (B.6) in the form 

Substituting Eq. (B.13) into Eq. (B.6), multiplying by W® and integrating over range Os^sl 
yields 


d 4 ^ 


\ 


W(x) + 0P,^)i!E = (Q,V) 

dxr 


where the inner product is denoted as 


<v v vj - f 


(B-14) 


(B.15) 


Thus the fundamental natural frequency of the beam can be obtained as 

/ \ 


2 


d 4x ¥ 


d% 4 




/(V,V) = -^ 


(B.16) 


where 
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N = 504( -36 + 60a - 60c 0 + 180ac 0 - 60c x + 120ac 1 

- 50c 2 + 90ac 2 -54c 3 + 84ac 3 - 60c Q c 3 + 60ac fl c 3 

- 300^3 + 3000^3 - 20c 2 c 3 + 20ac 2 c 3 - 15c 3 2 + 15ac 3 2 ) 

D = 28 + 252a + 168c 0 + 840ac„+ 1260c„ 2 + 1260ac 0 2 ^ ^ 

+ 120c t + 720ac t + 840c Q c 1 + 1680ac 0 c t + 210c 2 + 630ac 2 
+ 90c 2 + 630ac 2 420c Q c 2 + 1260ac Q c 2 + 252c t c 2 + 1008ac 1 c 2 

+ 84c 2 2 + 420ac 2 2 + 70c 3 + 560ac 3 + 252c 0 c 3 + 1008ac 0 c 3 
+ 1680^3 + 84000^3 + 120c 2 c 3 + 720ac 2 c 3 + 45c 3 2 + 315ac 3 2 

For the uniform beam discussed in section 5.1, which is clamped at one end and 
elastically supported at the other, the corresponding control parameters are y 2 =0, y 2 =0, kj=l, & 2 = 1 > 
Y 3 =l, y^= 1. If we substitute the middle point value of the identified boundary stiffness interval, 
i.e. jfc 3 =1221.65 and k 4 = 4.971, into Eq. (B.ll), the coefficients c, can be determined through Eq. 
(B.10) . Therefore, through Eq. (B.16), the nondimensional frequency by weighted residue 
method is obtained as \= 18.0637, which is associated with an error of 0.3539%. For the C-ES 
nonuniform beam, consider a special case a=0.5. For fcj=640.5 and k 4 =2A\5, Eq. (B.ll) results 
in a set of coefficients c, . Finally the weighted residues results in the following nondimensional 
frequency: X 2 =18.2561, which differs by 1.4228% from ^—18. 
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Fig. 3 The N-element beam which is elastically supported at both ends 
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Abstract 

In this study, the following two topics are considered for multi-span beams of both finite and infinite lengths 
with rigid transversal constraint and elastic rotational constraint at each support: (a) free vibration and the associated 
frequencies and mode shapes*, (b) forced vibration under a converted harmonic loading. The concept of wave 
propagation in periodic structures of Brillouin is utilized to investigate the wave motion at periodic supports of a 
multi-span beam. A dispersion equation and its asymptotic form is obtained to determine the natural frequencies. 
For the special case of zero rotational spring stiffness, an explicit asymptotic expression for the natural frequency 
is also given. New expressions for the mode shapes are obtained in the complex form for multi-span beams of both 
finite and infinite lengths. The orthoganolity conditions of the mode shapes for the two cases are formulated. The 
exact responses of both finite and infinite span beams under a converted harmonic loading are obtained. Thus, the 
position and the value of each peak in the harmonic response function can be determined precisely, as well as the 
occurrence of the so-called coincidence phenomenon, when the response is greatly enhanced. 

Introduction 

The model of a periodic multi-span beam with elastic supports is often utilized in 
engineering. For example, such a model is a reasonable approximation for a plate-like structure 
with parallel, regularly spaced stiffeners. The elastic supports may provide both the rotational 
and transversal restraints to the beam. Krein (1933) and Miles (1956) studied independently an 
N-span beam by using a finite difference approach, and established that the natural frequencies 
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fell into distinct bands with the same number of natural frequencies in each band as the number 
of spans. Lin (1962) generalized the finite difference approach to a multi-span beams with elastic 
supports. Abramovich and Elishakoff (1987) generalized Krein’s and Miles’s analyses to multi- 
span Timoshenko beams, taking into account shear deformation and rotary inertia. It was shown 
that the use of a finite difference approach might lead to computational difficulties and to 
inaccuracy in the determination of mode shapes of the system. 

Lin and McDaniel (1969) also used a transfer matrix formulation which is more 
convenient for the imposition of constraints at the supports. However, numerical difficulty may 
still arise when the number of periodic units in a structure is large . To overcome this difficulty, 
Yong and Lin (1989), and Cai and Lin (1991) transformed the state vector of displacements and 
forces into a vector of incoming and outgoing waves, and correspondingly transformed the 
transfer matrix into the wave-scattering matrix. By so doing, the computational efficiency and 
accuracy are greatly improved, especially when obtaining the dynamic response due to point 
excitation, because the calculation can be channeled in the direction of wave propagation. 

Mead (1970) made use of the concept of wave propagation in periodic structures 
originally due to Brillouin (1953) to analyze the free vibration of a multi-span beam of infinite 
length. Sen Gupta (1970) extended the analysis to finite multi-span beams and plates on rigid 
supports. In these studies, the wave propagation band and non-propagation band were studied 
in much detail. Sen Gupta (1970) also proposed a graphic method to determine the natural 
frequencies of the multi-span beams with rigid supports. In the framework of wave propagation, 
non-harmonic waves have to be decomposed into an infinite number of harmonic components 
in order to carry out the analysis of the forced vibration. This approach was used by Mead 
(1971) and Lin (1977) to obtain the response of an infinitely long multi-span beam to harmonic 
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excitation, as well as random boundary-layer pressure fields. In the actual calculation the infinite 
sum has to be truncated, and a large system of linear equations have to be solved to determine 
the unknown coefficients. 

It should be noted that in the forced vibration analysis of a periodic multi-span beam of 
finite length, multiple peaks occur in each wave propagation band. The number of peaks in each 
band is equal to the number of the spans. The computational effort becomes excessive when the 
number of span is large. Furthermore, an accurate position of each peak and its value are 
difficult to obtain. 

In order to circumvent the above difficulty, new expressions are proposed for the mode 
shapes of a periodic multi-span beam, based on the wave propagation concept, which can then 
be used in the forced vibration analysis. Since the transverse displacements within each span of 
the beam is related uniquely to the displacements at the two ends of the span, we may focus our 
attention only on the waves which propagate through each periodic support. Once these waves 
are determined, the motion of the beam between two neighboring supports can be computed if 
so desired. The dispersion equation which establishes the relationship between wave constant 
and frequency parameter is derived accordingly. The frequency parameters, wave constants and 
associated mode shapes for beams of both finite and infinite length can then be determined. The 
exact response of a multi-span beam to a convected loading is obtained for both the cases of 
finite total length and infinite total length. Furthermore, the locations of response peaks and their 
values can be precisely calculated, and the condition for the so-called coincident phenomenon can 
be predicted in exact terms. 


Free Vibration Analysis 
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Basic equations 


Consider an jV-span beam with uniformly spaced supports. It is convenient to write the 
equation of motion in terms of the local nondimensional coordinate ^ as follows: 

£/Wp (4) (i,0 + p AZ. 4 = 0 , ( P = 1,2 ,..., N ) (!) 

where is the transverse displacement in the 0-th span, and the local coordinate | is defined 


S = x/Z-(P- 1) , (P-1 )L*x*&L, 0^*1 W 

in which x is global coordinate and L is the individual span length. Assuming that the motion 
is harmonic 

I) * %(?)*-' , ( P - 1.2 N ) (3) 

where Wp(^) is the mode shape function associated with the p-th span, Eq.(l) can be reduced to 

w{ 4 \%) - = o W 


where 


^ pA c o 2 L 4 

{—El- 
is a nondimensional frequency parameter, and to is the sought angular frequency. 

It is assumed that each interior support provides a rigid constraint against transverse 
motion, as well as an elastic constraint against rotation, with a spring constant k (see Fig. 1). 
Thus the continuity conditions at each interior support are as follows: 
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( 6 ) 


%(1) = % +1 ( 0) = 0 
wtfi) = Wp^O) 

v^ p (l) = W^O) - Wp'(l) , ( P = 1,2 ,...,jV-1 ) 

where v=kL/EI. The first two conditions in Eq.(6) represent, respectively, the continuity of 
vertical and angular displacements. The last condition in Eq.(6) is the requirement of moment 
equilibrium at each interior support. The conditions at the two end-supports for a multi-span 
beam of finite length will be specified later. 

The mode shape that satisfies the first condition in Eq.(6) can be written as 

%(?) - A,ASA) * VU-SA) . ( P = 1,2, ...,N ) (7) 

where A p and 5 p are unknowns. Only the ratio Ap/B p is of interest in the free vibration case. The 
function XIA) is defined as 

ASA) - sin(X!) - Ji^sinhAS) ■ (8) 

sinh(X) 


Harmonic waves and the associated wave constants 

A simple wave of spatial sinusoidal variation cannot propagate along a multi-span beam, 
due to reflection at each support, giving rise to hyperbolic terms in the expression for the 
displacement. However, the concept of wave propagation can still be applied in the case of 
periodically supported beam, by focusing our attention on the waves which propagate through 
each support. The motion of a beam segment between two consecutive supports can then be 
determined from those of the two supports, if so desired. Since all the supports are assumed to 
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be transversely rigid, only the angular displacement at each support needs to be considered. Let 
0 p (t) be the angular displacement at the p-th support, and let it be represented in the form of a 
harmonic wave propagating through the p-th support, i.e 

0p(O =C e*-'-*) =©/'“' 

0 p = C e-** , ( P =0,1 ) 

where the nondimensional parameter p is known as wave constant, and is the amplitude of 
the propagating wave associated with the wave constant p. A positive p corresponds to a wave 
propagating in the positive x-direction, whereas a negative p corresponds to one propagating in 
the negative x-direction. The angular displacement function ©p is related to the mode shape 
function Wp(S=) as follows 

<«> I,.. - e s-l £ (10) 

w, '(?)!,., * V , (f> = 1,2,..., JV) 

obtained from the second condition in Eq.(6). 

The ratio A p /5 p in Eq.(7) will be determined for two special cases: The first case is 
associated with those mode shapes which are either symmetric or anti-symmetric with respect to 
the mid-point of each span. In such a case, the angular displacements at the two ends of a span 
are related as 

0 p-i = ( _1 ) f0 p » 5 = inte 8 er of 

where s=odd and s=even correspond to the symmetric and anti -symmetric mode shapes, 
respectively. In view of Eq.(9), the value of the wave constant for this case must be fi—ntK, 
implying a non-propagating wave or standing wave. The ratio A^JB^ is obtained by substituting 
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Eqs.(7) and (11) into Eq.(10) to yield 


= (-1 y +1 . ( 12 ) 

The second case is associated with those mode shapes which are neither symmetric nor 
anti-symmetric with respect to the mid-point of each span. Therefore, the wave constant fx is not 
an integer multiple of n, i.e. fx*rrm. This implies that there is indeed wave propagation through 
each periodic support of multi-span beam. By applying Eqs.(7), (9) and (10), we find that and 
Bp are given by 

A t -ii-mxje,., -/'(i,x)e,] 

, (13) 

b „ -il-rtUie,., */'(o.i)9 ( ] 

a - [/'(lA)P - [/'(o.MF - o . 

We note in passing that the first case corresponds -precisely to A = 0, associated with the same 
mode shape Wp(|) as that of a single-span beam with either two simply supported or two fully 
clamped ends. Return now to the second case, and substitute Eq.(9) into Eq.(13) to obtain 

A e = , n (14) 

Bp v\ m e~ ifl 

where 

n . (15) 

It is noted that the ratio A p /B p is independent of the span number 0 for both cases. This implies 
that one can choose a mode shape from any span as a reference, then the mode shape for the next 
span can be obtained by a phase shift. Thus, a general expression for the mode shape of a multi- 
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span beams may be written as follows 


wyi,/a) = (16) 

= [aj{%,\) * bj{ l-i,X)]e-> ( P- l) , ( P = 1,2, ...,N ) 

which is dependent on the span number P, the local coordinate |, the wave constant p and 
frequency parameter X. Here, a and b are obtained from Eqs.(12) and (14) 


a 



fi * mn 
= mn 



-e~ ifl r\* , fi * mn 
(-1 ) ,+l , fx = mn 


(17) 


Coefficients a and b are generally complex, while function /(') is real - Moreover, the span 
number p in Eq.(16) appears only in the exponential function. It will be shown later that this 
characteristic of mode shape is very useful, and it will be applied in the analysis of forced 
vibrations. 

Substituting Eq.(16) into the last condition in Eq.(6), the bending moment equilibrium, 
we obtain a dispersion relationship between fi and X as follows: 

cos 00 = F(X) , (18) 


where 

_ /'( 1,X) + / /2 (1,X) -/ /2 (0,X) (19) 

1 f'{ 0,X) 2f'(l,X)/'(0,X) 

Eq.(18) shows that the values of // and X must satisfy a certain relationship for the wave 
propagation. 

To examine the physical meaning of the dispersion equation, function F(k ) is plotted in 
Figs. 2(a) and 2(b). It is seen that F(k) has an oscillatory character; thus, each fi value 
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corresponds to multiple values of X. For the F(X) values between +1 and -1, the corresponding 
wave constants pi are real. This implies that there exists a non-zero phase difference between the 
motion in adjacent spans, and that the wave is propagating and the wave energy is being 
transferred from span to span without decay. The associated frequencies are grouped in 
distinctive bands, called the propagation bands. On the other hand, if the absolute values of F(X) 
are greater than 1, then pi is purely imaginary, indicating an exponential decay of wave motion 
from span to span. The corresponding frequencies are also grouped in distinctive bands, called 
the non-propagation bands. 

As shown in Figs. 2(a) and 2(b), the wave constant pi corresponding to the bounding 
frequencies of a propagation band must be an integer multiples of ji. At such a frequency , the 
motion of a multi-span beam reduces to a standing wave the same as that of a single-span beam 
with symmetric boundary conditions at the ends. The lower bounding frequency of the s-th 
propagation band is the same as the 5 -th natural frequency of a single-span beam with elastic 
rotational springs at the ends, whereas the upper bounding frequency coincides with the 5-th 
natural frequency of a single span with fully clamped ends. 

It should be noted that if pi is replaced by 

pi m = pi + 2 nm , ( m = ±1,±2,... ) (20) 

the dispersion equation, Eq.(18), remains unchanged. Thus the state of vibration of the system 
corresponding to a wave constant pi will be identical to the state corresponding to the other wave 
constant, namely pi+2nm. Therefore, if we want to have the one-to-one correspondence between 
the state of vibration of a system and the wave constant pi, the latter must be confined to a range 
of values of width 2ji. The range of pi values satisfying 
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( 21 ) 


(m- 1)ji < fi * mn 
-mi < fi s -(-m-l)jt , ( m = 1,2,3,... ) 

is known as the m-th Brillouin (1953) zone. For free vibration analysis, we may restrict to the 
first Brillouin zone ( m= 1 ) without loss of generality, i.e. 


0 < fi ss n , 
-n< fi £ 0 . 


( 22 ) 


We reiterate that a positive fi corresponds to wave propagation in the positive x-direction and a 
negative fi corresponds to one in the negative x-direction. 


Asymptotic dispersion relations and natural frequencies 

As seen in Eq.(18), the frequency parameter X is a multi-valued function of p. Let \{p) 
denote the X value in the 5 -th propagation band. For a large value of X», the following 
asymptotic approximation is sufficiently accurate 


FIX) - ( 1 + JL-)cos(X) - sin(X) , for Ui 
2X 


(23) 


which is obtained from Eq.(19) by letting tanh(x)-l and sinh l (x)-cosh l (x)-0. This 
approximation is compatible with the dynamic edge effect method due to Bolotin(1961) and 
Elishakoff (1976), and is remarkably accurate as shown in Figs. 2(a) and 2(b). Moreover, as the 
rotational spring stiffness v increases, the position of the lower bounding frequency of each 
propagation band moves toward the upper bounding frequency which is fixed. This implies that 
the multi-span beam structure becomes more rigid with larger v, as expected. In the case of v-0, 
the explicit asymptotic expression for the natural frequencies are obtained by combining Eq.(18) 
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and Eq.(23) 


, ( s = 1,2,.. .,00 ) (24) 

where s denotes the serial number of the propagation band. Thus, the natural frequencies of a 
multi-span beam can be determined readily from a given wave constant p which depends on the 
exterior boundary conditions of the entire system. 


\(ju) - (s - — )n + cos' 1 
4 


(_iy cosQf) 


The mode shapes of a multi-span beam 

It should be recalled that only the boundary conditions at the interior supports were used 
in obtaining an expression for W^,p,\). This implies that the mode shape given in Eq.(16) is 
valid only for a multi-span beam of infinite total length. For a finitely long multi-span beam, 
wave reflections occur at two exterior boundaries. Therefore, wave propagating in both positive 
and negative directions should be included in the analysis. The total angular displacement at the 
P-th support is now given by 


© p (a0 = + ©,(-a0 

= C e-w + , ( P = 1,2 


(25) 


where the positive and negative subscripts denote the two directions of wave propagation. Hence, 
the associated mode shape for a finite multi-span beam becomes 


W^,p,X) = C' + Ci Wpd.-^X) (26) 

where the p value will be chosen in the first Brillouin zone as defined by Eq.(22), and chosen 
to be positive without loss of generality. 

For an infinitely long multi-span beam, the wave constant p varies continuously over the 
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entire zone defined by Eq.(22). The associated frequency parameter X also varies continuously 
over the entire propagation band. For a finitely long multi-span beam, however, the wave 
constant fx and the associated frequency parameter X take on discrete values. The number of the 
discrete values for n or X in each propagation band is the same as the number of spans. These 
discrete wave constants are determined by imposing the boundary conditions at the exterior ends 
of the entire beam. Referring to Fig. 1, the boundary conditions at the exterior ends of the beam 
are 


vX(0„w,X) - <(0,/i,X) , 


K L 


v^(l,^,X) = - Wjf , v 0 = -2_ 


EI 


v * = 


k N L 

El 


(27) 


where v 0 and \ N are the non-dimensionalized rotational spring constants at the left and right ends 
of the N multi-span beam, respectively. A vanishing v corresponds to a simple support, and an 
infinite v to a clamped support. 

The mode shape of a finitely long multi-span beam can be rewritten in abbreviation as 

follows 


W'JI) - = s, JjX) ♦ v, /;(!-!) . 


(28) 


where 


(29) 


a R . = A. [a. s + p.a 

P.; / L i > > J ’ 

v = . 

P.y y L i i > J 

and where the subscript (3 denotes the (3-th span, fx-. is the wave constant corresponding to X ; , and 
A ; and T ; are the unknown constants to be determined by imposing the boundary conditions at 
the exterior ends. 
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Examples 


The wave constant p , frequency parameter X and mode shapes for an Af-span beam will 
be evaluated in detail for the following three cases. 

(a) Case v,=v N =v/ 2 

In this particular case the rotational spring constants at both ends of the multi-span beam 
are equal to one-half of those at the interior supports. The boundary conditions at the exterior 
ends are 


lw[(0,p,X) = w"(0,p,X) , 
lw'(l,p,X) = X) . 


(30) 


Using Eqs.(18), (26) and (30), we obtain, after some algebra, an equation for p as follows: 

sin(uAf) = 0 • 

The possible values of fx in the first Brillouin zone are 

p. = Ln, (j = 0,1,2,...,AT) . ( 32 ) 

' N 

As seen in Fig. 2, the values p= 0 and //= n are associated with the bounding frequencies of the 
propagation bands. In the case of an odd-numbered propagation band, p is equal to zero at the 
upper bound and to k at the lower bound. The opposite is true for an even-numbered 
propagation band. Moreover, the lower and upper bound frequencies are the same as a single- 
span beam with elastic supports of rotational spring constant of value v/2, and with fully clamped 
supports, respectively. To incorporate the above features, Eq.(32) for p is modified to read 
where the subscripts s and r denote, respectively, the s-th propagation band and the r-th 
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/V(,-1)*T | 2^ ^ ^ N j 31 ’ (33) 

( 5 = 1,2,.. .,» , r = 1,2 ) 

frequency within each band. Then the frequency parameters X ; _ (j . 1 y V+r is numbered in an 
increasing order of j. 

The mode shape W^/l) should be taken as the real part of Wp(|,/Zy,X ; ), and the 
coefficients A and T ; in Eq.(29) should assume the values of 



(34) 


( b ) Case v # =v N =°o 

In this case, the boundary conditions at the extreme ends are 

= ® 0 (fi)L = 0 , 
wft 1,A*,X) = = 0 • 

Eq.(31) remains valid; however, the serialized version now reads 


- (-I ) 1 ^ 


IJt , 


( s = 1,2,...,® , r — 1,2,... JV ) . 


(35) 


(36) 


The wave constant pc= 0 and fi=n correspond to the upper bound of an odd-numbered and an even- 
numbered propagation bands, respectively, contrary to case (a). In this case the mode shape takes 
the imaginary part of Wp(|^-,Xy) and the coefficients in Eq.(29) are found to be 
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1/2 , 

t*i = kn 

i/2 , 

fi. * kn , 

1 , 

It 

5C 

-1 , 

/a. * kn . 


( 37 ) 


(c) Case v,=v/2 and v^oo 

In this case, the left end of the multi-span beam is constrained by a rotational spring of 
stiffness constant v/2, while the right end is clamped. The corresponding boundary conditions 
read 


0,//,X) = W"(0,^X) , 
= 0^)1 = 0 . 

Analogous to cases (a) and (b), the following equation for fi is obtained 


from which 


cos (fiN) = 0 , 


(38) 


(39) 


*(-!)* 

( 5 = 1,2,.. , r = 1,2 ) . 

Note that neither zero nor n is a solution of the above equation; thus a standing wave does not 
exist. The mode shape for this case is described by the real part of W p (|,yU ; -,Xy) with coefficients 
specified in Eq.(34). If the left end of the multi-span beam is treated as being clamped (v 0 =o°), 
and the right end is treated as being elastically constrained by a rotational spring stiffness of v/2, 
then the mode shape is described by the imaginary part of Wp(|,^,X ; ) with coefficients given in 
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Eq.(37); while Eqs.(39) and (40) remain unchanged. 

Tables I-III list the frequency parameters in the first two bands of a six-span beam 
evaluated by using both the exact and asymptotic formulas for the three sets of boundary 
conditions at the extreme ends. It can be seen that the exact and asymptotic solutions differ by 
less than 0.4 % in the first band, and they are almost identical in the higher bands (si2). 

The normal modes associated with the frequency parameters in the first band are 
illustrated in Figs. 3(a), 3(b) and 3(c), respectively. The solid and dash lines correspond to the 
cases of v=0 and v=10, respectively. Figs. 3(a) and 3(b) portray the mode shapes for the two 
sets of exterior supports, namely the set of Vq^Vj^v/2 and the set of ‘v 0 =v N = l ^>. In these two cases 
the mode shapes are either symmetric or anti-symmetric with respect to the mid-point of the 
multi-span beam due to symmetric boundary conditions at the two extreme ends. Fig. 3(c) 
illustrates the mode shapes for the case v 0 =v/2 and Vy=«>, and they are neither symmetric nor 
anti-symmetric, as expected. 


Forced Vibration Analysis 

In the preceding section, the frequency parameter X, wave constant p. and the associated 
mode shape have been determined for a multi-span beam of finite or infinite total length. In this 
section, the exact analytic harmonic response of such a beam subjected to a convected harmonic 
loading is obtained using the normal mode approach. Furthermore, both the location and the 
magnitude of the peak response can be determined in advance. Thus, the important coincident 
phenomenon can be investigated in exact terms. 
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Orthogonality conditions of mode shapes 
(a). Multi-span beam of finite total length 


Consider two normal modes of a multi-span beam satisfying the following equations 

"ft’d) -K (41) 
lift’d) - K liftd) - o . 

where the first subscript (3 denotes the fi-th span, and the second subscript, j or k, corresponds 
to the serial number of a natural frequency. Multiplying the first equation in Eq.(41) by t (|) 
and the second by Wp ; (i), and integrating the difference between the two resulting expressions 
over the total length NL of the JV-span beam, we obtain 



[liftd) lift’d) - »ftd) lift’d)]^) 






( 42 ) 


— — 

Integrating the product W p t (l)Wy (|) by parts 




(43) 

1 

-iiftd)>iftd)]!, + f "ftd) lift’d) • 


Eqs.(42) and (43) can be combined to yield 


s 


E [i*ftd)ii£d) 


nftdJwftd) * wftdjwftd) - %'»(?)%.,<?)]» 

(44) 

> d‘ - K)f. fiiftd)iiftd)‘d . 

f-i 4 
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Finally, taking into account the continuity conditions Eq.(6) at each interior support, we obtain 


- KkWKiW * 

-[^,,(0)^7°) - 


* »£,(0)tf w (0) - ffjoiffJO)] 



(45) 


It may be noted that the left-hand side of Eq.(45) vanishes for any set of homogeneous boundary 
conditions of the form 


aW(x) + bw"(x) = 0 


(46) 


or 


cW'(x) + dW\x) = 0 ^ 47) 

at the two ends, where a, b, c, d are constants. The idealized boundary conditions, such as 
clamped-free, simply-simply supports, and so on, are special cases. Eq.(46) corresponds to a 
transverse elastic support, and Eq.(47) to a rotational elastic support. Thus, we obtain an 
orthogonality condition for normal modes of an AT-span beam as follows 


p-i i 


(48) 


where 6 ;Jt denotes the Kronecker delta and y/ is defined as follows: 

p-i { 


(49) 


Note that the span serial number p and the local coordinate are separable in the 
expression for the mode shape given in Eqs.(28) and (29). Indeed, p appears only in the 
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exponential functions not involving Therefore, integration over the entire length of a multi- 
span beam can be carried out with respect to the local coordinate and then summed over the 
span serial number p. It will be shown later that these properties can be used to advantage in 
reducing the computational efforts when evaluating the dynamic response of the system. 

Eq.(49) may be rewritten as follows: 


Yy - E - CA(*,) 

P-i i 



(50) 


where I k (-) and / 2 (-) are integrals defined as 




-ill * J-SUX2X) - sin ^> . 

2 1 2X sinh 2 (X) 


1 + sinh(2X) 

2 X 


(51) 


w • 


1 

|Ai)Ai 




= J_ I -cos(X) - -Lsin(X) + s * n i M 
2 1 X sinh 2 (X) 


cosh(X) + — sinh(X) 
X 


(52) 


and where C “ and C~ can be obtained from the following more general expressions 


JV 


(53) 


Cjk = 52 ( a p./ a p.* + ^p.j'^p.* ) 

p - 1 

= A.A t {(a.a t + + r.r t (a/a; + b/b t ’) S*(^. + {i k ) 

+ r. ( af a k + by* b k ) S w *(/i. -/* t ) + ( a- a** + b b t * ) ~M k ) } , 

In Eqs.(53) and (54), an asterisk denotes the complex conjugate. T’s and A’s are defined by 
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+ fl eA;) 


( 54 ) 


c i ■ E ( sA* 

P-1 

- A ; A, {(«,*, * a„b^S^t <,) * TS t (a-b; * at b')S^^) 

* r „ ( a , b i * a k b i') S ^t i i~f‘t) * r , ( b t * a t b j )S^(/i i -/t l ) } 

i, %<S) -Jfe{w f (?)} 

- 1 , W t (%) «/m{W ( ($)} . 

and function is given by 


A = 


1/2, W p (|) =to{W p (i)} 

-in , W{%) =/m{W p (i)} 


r = 


S M = 


N 


E 


e -*>(P- 1) 


1 - <? 
1-e'* 

N 


H * 2mn 
fi = 2mn . 


(56) 


(b). Multi-span beam of infinite length 

In the case of a periodic beam of infinite length, the orthogonality condition of normal 
modes can also be derived by using a similar procedure. However, it is no longer necessary to 
impose any boundary conditions at the exterior supports. The mode shape Wpf^^UjXXu)] given 
in Eq.(16) is now applicable throughout the entire length. Eqs.(41) through (44) still hold, except 
that the finite sum is replaced by an infinite sum. By taking into account the continuity 
conditions, Eq.(6), at the interior supports, it is easy to show that the left-hand-side of Eq.(42) 
vanishes, i.e 
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2: #y){(e^ -e-'O^[0,/i,XW]^[l/W)] (5?) 

+ («-* - e^)W l '[l,Ai,X/^)]W r 1 "[0y,X / (/i')] } - 0 , 
where 6(-) is Dirac’s delta function, and where use has been made of the identity 

oo 

£ e -,XP-i) = 2* 60*) . ( 58 ) 

P m -00 

The orthogonality condition of the mode shapes for an infinitely long multi-span beam read 

« i 

£ fW’ p [|,A*.W]W r p[l.A*'.W)]^ = Yo[i“A»]6(^V)^ , ( 59 > 

P— « * 

where 6^ is the Kronecker delta, and y 0 2 [-] is defined as follows: 


= ( |aM*l 2 ) /,[*•>) 1 *2 Re{ab'} /,[X,M] 


(60) 


in which integrals If-) and / 2 (-) are given in Eqs.(51) and (52), respectively, and a, b are defined 
in Eq.(17). 


Responses of multi- span beam under the connected loading 

Let us consider the forced vibration of a multi-span beam with damping. The equation 
of motion in the local coordinate system is given by 

£/L- 4 y p (4) (i,0 + c^(l,0 + pAy p (i,0 • p f (hQ , (61) 

( p = ; 0 s | * 1 ) 

where c = damping coefficient, p^t)= transverse pressure per unit length, and N=l and N=M 
for a finitely long beam, whereas N_=-°° and N=<p for an infinitely beam. Assuming that the 
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excitation and response are harmonic in time, we have 


PW* ( 62 ) 

Function Fp(|) will be referred as the harmonic- response function in what follows. For a 
harmonic loading convected over the beam at a velocity co Llft f 

P p (i) = P 0 e ‘ ifl ' (i * , ( p = 1,2 JST ) ( 63 > 

where P 0 is the amplitude, and is the wave constant of the loading. 


(a) Multi-span beam of finite total length 

First, let us consider the case of an jV-span beam. We expand Yp(|) and P p (%) in terms 
of the normal modes of the system as follows 


y,<X) - E . 

y-i 

00 

f = 1,2,..., AT ) 
y-i 


(64) 


The relationship between coefficients c ; and d) can be found by substituting Eqs.(62) and (64) into 
Eq.(61) and using Eq.(41) to obtain 


E ( - K - = 


El 


E W$> 


( 65 ) 


7-1 


where X 4 =pAL 4 (a 2 /EI is a nondimensional loading frequency parameter, and Z=cL 2 /(pAEIf 1 is a 
nondimensional damping parameter. Comparison of coefficients on the two sides of Eq.(65) 
yields 
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Cj - , 

where is the frequency response function given by 


( 66 ) 


- V * • 

Hence, the harmonic response function is obtained in the following form 


(67) 


iys> - e w , c p - 1,2 jv ). 


( 68 ) 


The coefficients <i ; are obtained by applying the orthogonality condition for the mode shapes, 
namely Eq.(48), to the second equation in Eq.(64) to yield 


where 


d. = 

j 


i h 1 _ 

4E K/i)p p (i)^ 

Y j o- 1 i 


_ Po 


= -JL[ E - + E/e-^gXV^)] , 


Y; 


(69) 


V = 


2 + 

in which sgn(-) denotes the sign function. 


r j a ’ s ^ f -^) ] , 

(70) 

T iK s M~p) 1 > 

(71) 

= U - V , 

(72) 

sinh(X) 

+ ifi f ) e x ~‘* r - 2\] , 

(74) 


3-23 



r [(W,)^'* 5 + ( X.-A)e“' (k ^ ) - 2X] , |^|*X 

. -i 7\ 1 * J 


2(rf-V) 


(73) 


--isgn^) * __i_^ { 1 - } , |*|-X 


(b) Multi-span beam of infinite length 

The dynamic response of an infinitely long multi-span beam subjected to a convected 
harmonic loading can also be evaluated in a similar way. Let us expand both the harmonic 
response and the loading function in the mode shapes of such a beam 



where s denotes the serial number of a propagation band. Multiplying both sides of the second 
equation in Eq.(75) by VK p *[^, J u',X J ,(//)] and performing integration over the length of the entire 
beam, we obtain upon applying the orthogonality conditions, Eq.(59), 


= — r i i; | . (76) 

2 Jt Yo| |>,\( a 0 ] J 


For an external excitation in the form of Eq.(63), the numerator in Eq.(76) may be simplified to 
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(77) 


oo 1 

E K(S) W fT[S,M>)]^ = 2nb( Mr riP 0 D[fi,\(n),ti f ] , 

P*-® TJ 

D[n,\(n),n f ] = a m [p,X s (fi)] g[\(/i),fi f ] + b’[fi,\(fi)]e~ ,flr , 

in which a(-) and b{-) are defined in Eq.(17), and g(-) is given by Eq.(72). Eq.(66) is still valid 
for cjjx) and djji). Hence, c s (/i) may be expressed as follows 

ejji) - mhiri.M,] ■ < 78 ) 

■ifrA/M] 


The harmonic response function for the infinitely long multi-span beam can be obtained by 
substituting Eq.(78) into the first equation in Eq.(75) to yield 


'- 1 4 * Yo [/*,\(/*)] ^ 

( p = 0,±l,±2,...,±oo ) . 

Here, the integration range R must be chosen from the particular Brillouin zone defined in 
Eq.(21) which includes the loading wave constant fip This is always possible since union of all 
Brillouin zones constitutes the entire one-dimensional space. Note the presence of a Dirac’s delta 
function in Eq.(79) with an argument It implies that only a group of propagating waves 

associated with /^contributes to the response. Carrying out the integration in Eq.(79), we obtain 


- E 





( 80 ) 


( p =0,±l,±2,...,±oc ) 


where y 0 \-), D(-), H(-) and W p (-) are given by Eqs.(60), (77), (67) and (16), respectively. 
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Results and discussion 


For a multi-span beam under harmonic excitation, large response is likely to occur if the 
excitation frequency is within a propagation frequency band. When the excitation is convected 
along the beam, the response can be further amplified due to the so-called coincidence effect. 
For an JV-span beam the magnitude of the harmonic response function may have as many as N 
peaks in each propagation band. The possible location of a peak in a propagation band for |ip(^)| 
can be determined from the condition 




1 , 2 ,... ) , 


(81) 


when the denominator in expression (67) for |ff(X^X^)| has a minimum magnitude, whereas the 
magnitude of harmonic response function lp(^) is dominated by the term with |/7(X ; ,X / )|. For an 
infinitely long multi-span beam, however, there is only one group of propagating waves, whose 
wave constant coincides with that contributes to the response, as it can be seen from Eq.(79). 
Therefore, there is only one peak at 



( 5 = 1 , 2 ,... ) 


(82) 


appearing in each propagation band. 

Fig. 4(a) portrays the harmonic response at the mid-point of the second span of a four- 
span beam. It shows that there are four peaks in the first propagation band. In contrast, there 
is only one peak in the first band for the infinite multi-span beam shown in Fig. 4(b), and the 
response is considerably magnified due to the coincidence effect. 

The effects of damping on the harmonic responses are also shown in Fig. 5(a) and 5(b) 
for a four-span beam and an infinitely long multi-span beam, respectively. It can be seen that 
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the values of the peaks in each propagation band is reduced with a larger damping coefficient. 
Hence, the profile of each peak becomes flatter as damping increases. 

The results obtained from the present approach and Mead’s approach (1971) are compared 
in Figs. 6(a) and 6(b) for an infinitely long beam with an evenly spaced hinge supports (v=0). 
The dash line represents a 20-term approximation in Mead’s formulation, whereas the solid curve 
represents a one-term approximation by the present approach in the first propagation band. The 
results are seen to be very close. The two results become indistinguishable in the first three 
propagation bands, when a twenty-term was used in the present approach, as shown in Fig.6(b). 
Most significantly, the present approach has the advantage in that the location of each possible 
peak in each propagation band can be determined; thus, the value of each peak can be evaluated 
precisely from Eq.(80). 


Conclusions 

The free and forced vibrations of periodically supported multi-span beams are studied in 
the this paper. Both the case of finite total length and the case of infinite total length are 
considered. The wave propagation concept is applied in the analysis of free vibration of the 
beam systems. The dispersion equation and its asymptotic form are derived from which the 
natural frequencies can be determined for a given wave constant. An explicit asymptotic 
expression for the natural frequencies is also proposed for the specific case of zero rotational 
spring stiffness. It is shown that the agreement between the asymptotic and the exact frequencies 
is excellent. The mode shapes of free vibration are obtained in the complex form. In these mode 
shapes the span serial number and the local spatial coordinate are separable; thus, an integration 
over the entire length is reduced to one within a single span, and a summation over the serial 
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numbers of the spans. It is shown that the use of these mode shape expressions can greatly 
reduce the computational efforts in the forced vibration analysis. 
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Table II. Frequency Parameters of Six-Span Beams with 
Rotational Spring Parameter v ( Case v 0 =^*,=00 ) 


— 

w-kLIEl 

0 

2 

200 



Exact 

Asymp. 

Exact 

Asymp. 

Exact 

Asymp. 



Eq.(18) 

Eq.(24) 

Eq.(18) 

Eq.(23) 

Eq.(18) 

Eq.(23) 


Frequencies in 

3.261 

3.267 

3.491 

3.491 

4.647 

4.630 


the first band 

3.556 

3.566 

3.729 

3.730 

4.663 

4.646 




3.927 

3.927 

4.042 

4.037 

4.685 

4.668 



4.298 

4.288 

4.362 

4.351 

4.724 

4.690 




4.601 

4.586 

4.623 

4.607 

4.724 

4.706 

— - 


4.730 

4.712 

4.730 

4.712 

4.730 

4.713 

— 


6.410 

6.410 

6.536 

6.536 

7.720 

7.720 


Frequencies in 

6.708 

6.707 

6.802 

6.802 

7.746 

7.746 


the second band 

7.069 

7.069 

7.134 

7.134 

7.781 

7.782 



7.430 

7.430 

7.468 

7.468 

7.817 

7.818 



7.727 

7.728 

7.740 

7.741 

7.843 

7.844 



7.853 

7.854 

7.853 

7.854 

7.853 

7.854 
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Table III. Frequency Parameters of Six-Span Beams with 


Rotational Spring Parameter v ( Case v 0 =v/2, ) 


v=kL/EI 

0 

2 

200 


Exact 

Asymp. 

Exact 

Asymp. 

Exact 

Asymp. 


Eq.(18) 

Eq.(24) 

Eq.(18) 

Eq.(23) 

Eq.(18) 

Eq.(23) 

Frequencies in 

3.173 

3.175 

3.422 

3.421 

4.643 

4.626 

the first band 

3.393 

3.403 

3.596 

3.598 

4.654 

4.637 


3.738 

3.743 

3.881 

3.879 

4.674 

4.656 


4.116 

4.111 

4.205 

4.196 

4.697 

4.679 


4.463 

4.451 

4.505 

4.491 

4.717 

4.699 


4.696 

4.679 

4.702 

4.685 

4.728 

4.711 


6.317 

6.317 

- 6.456 

6.456 

7.713 

7.713 

Frequencies in 

6.545 

6.545 

6.656 

6.656 

7.731 

7.732 

the second band 

6.885 

6.885 

6.964 

6.964 

7.763 

7.763 


7.252 

7.253 

7.304 

7.304 

7.800 

7.800 


7.592 

7.592 

7.617 

7.618 

7.832 

7.833 


7.820 

7.820 

7.823 

7.824 

7.851 

7.852 
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Fig. 1. Multi-span beams with elastic rotational spring at each support 
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Fig. 2. Dispersion relationships 
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v = 10 



3.897 


3.950 


4.094 


4.293 


4.501 


4.66 


3(a). v 0 =v/2 and v N =v/2 

Fig. 3. Normal modes in the first band of a six-span beam ( — for v=0; — for v=10 ) 
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v = 0 ; v = 10 



K = 3.173 ; 


X* = 3.393 ; 


X, = 3.738 ; 


X 4 = 4.116 ; 


X, = 4.463 ; 


X* = 4.696 ; 


3(c). v 0 =v/2 and v N =°° 


3.911 


4.013 


4.189 


4.399 


4.593 


4.713 
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IW)W/£7)| 


p=2 |=0.5 ^=4.14 


Fig. 4(a). Nondimensional harmonic response of a four-span beam vs 
loading frequency parameter X f ( — Using normal modes of 
the first band; — Using normal modes of the first two 


0.5 4.14 


Fig. 4(b). Nondimensional harmonic response of an infinitely long 
multi-span beams vs loading frequency parameter ( — 
Using normal modes of the first band; — Using normal 
modes of the first two bands ) 









|K p (£)/(PoL 4 /£/)| \Y0)t(P£ A /El)\ 



X/ji 


Fig. 5(a). Nondimensional harmonic response of a four-span beam vs 
loading frequency parameter with different damping 
coefficients ( fJ=2, £=0.5, /^=4.14 ) 



Xyfo 


Fig. 5(b). Nondimensional harmonic response of an infinite long 
multi-span beam vs loading frequency parameter X^ with 
different damping coefficients ( 0=2, £=0.5, ^=4.14 ) 
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p=2 |=0.5 ^=4.14 


Nondimensional harmonic response of an infinitely long 
multi-span beam in the first band ( by present approach; 
— by Mead’s approach ) 



Fig. 6(b) Nondimensional harmonic response of an infinitely long 
multi-span beam in the first three bands ( — by present 
approach; — by Mead’s approach ) 
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VIBRATION OF TWO-DIMENSIONAL GRILLAGE 


L. P. Zhu, Y. K. Lin and I. Elishakoff 
Center for Applied Stochastics Research and Department of Mechanical 
Engineering, Florida Atlantic University, Boca Raton, FL 33432 


Abstract 

In this paper, an exact solution for free vibration of a two dimensional periodic 
grillage structure is obtained by using wave propagation concept. Each periodic node (or 
support) of such a grillage structure is constrained by a rigid transverse support, and two 
elastic rotational springs in two orthogonal directions. The rotational springs in each 
direction are identical. The four boundary edges of the grillage as a whole are either 
simply supported or clamped. The wave motions at all periodic nodes are investigated 
and the associated dispersion equation relating the natural frequency and wave constant 
is derived. It is shown that there exist both bending and torsional motions in a grillage 
structure, which are coupled between members in two perpendicular directions. General 
expressions of the synchronous mode shapes and non-synchronous deflection shapes are 
obtained for both bending and torsional motions in two directions. These can then be 
applied in formulating forced vibrations of the grillage structure. An example of a 
grillage with four-spans in each direction is given for illustration. 
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1 . 


Introduction 


A grillage consists of two sets of mutually orthogonal beams rigidly connected at 
the intersections where additional elastic constraints may be imposed in rotational and/or 
transverse directions. The properties and spacing of the beams, and the constraining 
springs in each direction are spatially periodic; however, they may differ in two different 
directions. 

Grillages find wide applications in many forms of construction. A steel or 
concrete deck is often placed on a grillage to transfer the imposed loads. There are 
several methods for the analysis of a linear elastic grillage. The general method of 
statically indeterminate structures is applicable in the case of static loading, but it is 
extremely tedious except when the number of intersections is small. If the space between 
two neighboring beams is small then approximation of the overall grillage by an 
orthotropic plate is reasonable, for which the appropriate differential equation can be 
found (e.g. Timoshenko 1959). Ellington and McCallion (1957) proposed a more accurate 
method to obtain deflection under static loads by using a finite difference formulation, but 
they assumed that the beams did not resist torsion. 

There have been rather few papers on the vibration of grillages. In their 
pioneering work, Ellington and McCallion (1959) assumed that the mass of the beams 
was concentrated at the nodes (intersections) of the grillage. Again, by neglecting the 
torsional stiffness of the beams, and using a finite difference equation approach, they 
obtained a frequency equation for a grillage with various boundary conditions. Wah 
(1963) improved this method by considering the torsional stiffness of the beams with both 
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distributed and lumped masses. It was noted in his work that the derivation of the 
frequency equation required an admissible transverse displacement at each node which 
satisfies the boundary conditions. Unfortunately, such an expression is not easily 
determinable except for the case of a rectangular grillage simply supported on all four 
exterior edges. For other support conditions, only the case of zero torsional stiffness was 
discussed in his paper. Recently, Dinkevich (1993) considered a grillage with a 
transverse spring at each support. A frequency equation was obtained for the case of 
simply supported edges and negligible torsional stiffness. 

In this paper, both wave propagation in a periodic grillage of infinite size and free 
vibration of a grillage of finite size are investigated. Exact solutions for both cases are 
obtained by using the wave propagation concept. Every periodic node (or support) is 
constrained by a rigid transverse support and two elastic rotational springs in two 
orthogonal directions. The rotational springs are identical in the same direction, but they 
may be different in different directions. In the case of a grillage of finite size, the four 
boundary edges are assumed to be either simply supported or fully clamped. Both 
bending and torsional motions of each beam are taken into account in the analysis. The 
bending motions of the beams in one direction are coupled with the torsional motions of 
the beams in the other direction. The wave motions at all periodic nodes are investigated, 
and a dispersion equation that relates the natural frequency and the wave constant is 
derived. The general expressions of the mode shapes for both bending and torsional 
motions are obtained, which can be applied in the analysis of forced vibrations. An 
example of a grillage with four-spans in each direction is given for illustration. 
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2. Wave Propagation in an Infinite-Size Grillage 

2.1 Basic Equations 

Consider a two-dimensional grillage structure of an infinite size, composed of two 
sets of multi-span beams in two perpendicular directions. It is assumed that each support 
provides a rigid constraint against transverse motion, as well as elastic constraints against 
rotations about two axes, with elastic constants k x and ky , respectively, as shown in Figs. 
1(a) and 1(b). Although the transverse and torsional motions are uncoupled for a 
prismatic beam in one direction, they may be coupled in two sets of beams in the 
perpendicular directions. Let us focus our attention on node/support (a,f3) as shown in 
Fig. 1(a). The displacements in the four beams sharing the same support (a, (3) are 
denoted by the symbols as shown in Fig. 1(b). The equations for the bending motions 
of the beams in the x and y directions are, respectively, 


i 4 JSf fl (S,f) 4 d 2 X a a (|,0 

apV5 : + m l* ap ^: = 0 

d * dt 


( 1 ) 


D 




* m , L , 


, d’r.a.O 


~d? 


= 0 


( 2 ) 


where X a p ,Y^ = transverse displacements, m x ,m y - masses per unit length, L x , L y = span 
lengths and D x , D y = bending rigidities of the two beams in the x and y directions, 
respectively. 

The equations of the torsional motions are 
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( 3 ) 


|2 e a f »(£’0 


- n y L y 


2 d 2 QJK,t) _ 


= 0 


dV 


*Mi.o _ mT 2^ afl (bt) =ft 

d% 2 dt 


( 4 ) 


Here, and 0 are, respectively, the torsional displacements of the beams in the x and y 
directions, n x , n y = polar mass-moments of inertia per unit length, and C x ,C y = torsional 
rigidities. The effect of warping of a cross-section due to non-uniform torsion is 
neglected. The positive directions of and 0 are determined by the usual right-hand rule. 
The local coordinates t; and l, in the two directions are defined as follows: 


§ = x/L x - (a - 1) ; (a - 1 )L x s x s aL x ; Os^sl 
£ - y/L y - (P - 1) ; (P - 1 )L y sy*$L y ; O^sl 


For an infinitely large grillage, there are no exterior supports. Continuity 
conditions must be met at each support. For the transverse motion of the beam, they 
require that 


^1.0 1, .i -jw?.oi s .o- 0 

XLp (1,0 I,.! -*.'*1*5*0 U.o = ®.*0*, 


D 


*£ip(5-OI,.o 


+ c 




( 6 ) 


and 
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( 7 ) 


i t „ - n,.i«.o u - o 
yUm i 5 ., - n',.i(6.«) i t „ - 

o,[nV,(5>‘)i t .o - nVt-oit.i] 


- C 


^a-l^O |,.o - *U(i.O U-l] = * ^(0 


where a prime denotes differentiation with respect to % or For the torsional motion, 
the boundary conditions read 


— 

e aP (^0U 



= © ,(0 

( 8 ) 





- *„p(0 

( 9 ) 


in which ©<^(0 and v P a p(f) are, respectively, the rotations at the (a,P)th support about two 
horizontal axes as shown in Fig. 1(b). 

As seen from Eqs. (1) through (9), the ©(f) motion is related to the bending 
motion of the beams in the x direction and the torsional motion of the beams in the y 
direction, whereas, the ^(t) motion is related to the bending motion of the beams in the 
y direction and the torsional motion of the beams in the x direction. It is, therefore, 
possible that the vibration of a grillage structure can be decomposed into two uncoupled 
states in terms of the ©’$ and W’s, respectively. In what follows, we shall discuss the © 
motion only. The motion is similar. Assume that the motions is harmonic 
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( 10 ) 




0 a p(M =0^(0^“' 

0 o P (O = ©ap 6 '"" 

Then, Eqs. (1) and (3) become 

*Z (5) - *4* P (S) = o 

5^(0 ♦nie.^) =o 

where Xq and rj e are nondimensional frequency parameters given by 


f m L* 

1 

T 

/ r 2 \ 

n L v 

x x O) 2 

n = 

y y (U 2 

D 

> Me 

c 

x ) 


* / 


Solutions of Eqs. (13) and (14) are obtained as 

6« ( (?) * C =f ,sinr, 0 6 * D J( sinti s (l-:;) 
where function / x (-) is given by 


(ID 

( 12 ) 


(13) 

(14) 


(15) 


( 16 ) 

(17) 


A(i) = sinX| - 4^-sinhXJj (18) 

sinhX 

By applying the continuity condition at the node (a,P), the coefficients in Eqs. (16) and 
(17) are determined in terms of the angular displacements at the nodes (a-l,|3), (a,|3-l) 
and (a,|3) as follows: 
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A 


“P 


5 a P = 


A(^-e) 

L 

x 

AW 


KC 1 )®., -/4(°) 0 a-ip] 

[/4(°)®a P - <(!)©«- ip] 


( 19 ) 


'aP 


a P 


simi € 


sm»i e 


( 20 ) 


where A(-) is given by 


7 sinX„ 

A(X e ) = 2?4^_^_( 1 - cos*. coshX ) (21) 

sinhX e w 

It can be seen from Eq. (17) that sint] e = 0 corresponds to a single span beam with 
clamped supports for the torsional motion. 


2.2 Wave Motions and Dispersion Relation 

Similar to the case of one-dimensional multi-span beam discussed in previous 
paper [9], only the motions at the periodic supports (nodes) of the structure are to be 
investigated. Note that the motions at all periodic nodes (supports) of the grillage are not 
independent, but they may belong to that of a two-dimensional harmonic wave, which 
may propagate through all the nodes of the grillage. 

For the grillage structure presently considered, only rotational displacements are 
possible at the supports due to the transversely rigid constraint. Let 0 ap (r) be the angular 
displacement at node (a,P) shown in Fig. 1(b), which can be represented by a two 
dimensional harmonic wave propagating through node (or support) (a,P) as 
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0 .(0 -c 9 e‘'“ or 


© . - C e 

a p W 


( 22 ) 


in which and are the components of a wave constant vector in the x and y 
directions, and C e is the amplitude of the wave motion. 

If this two dimensional harmonic wave can propagate through the periodic nodes 
of the grillage, the frequency parameters Xq and tJq, which are functions of natural 
frequency to and wave constants must be related in a dispersion equation as 

follows 


2y x, 


sinhX 0 ( cosX 0 - cosp. Q x ) - sinX e ( coshX 0 - cosp e x ) 


* e ' 


cosX 0 coshX 0 - 1 


(23) 


2(l- 


Y l lie 


cosri e - cos^ 0 y 
single 


+ v =0 


where v x and y z are nondimensional parameters given by 


= k 


/ 

\ 

_ 1 


\ 

k 2 

C 2 

y 

T D 

X 

D 2 X 

Cy\ 

^ 1 

^1 



II 

L ’l 


1 

'7 


0 * Y, * 1 


(24) 


The dispersion equation is obtained by substituting Eqs. (16) through (20) into the third 
condition in Eq. (6). The first and second terms in Eq. (23) indicate, respectively, the 
contributions from the bending motion of the beams in the x direction and the torsional 
motion of the beams in the y direction, and the last term is attributable to the external 
rotational spring k x . When y x = 0 (or D x = 0), the problem reduces to a purely torsional 
vibration of a multi-span beam with periodic elastic rotational constraints in the y- 
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direction. On the other hand, when y, = 1 (or C y = 0), the problem becomes a purely 
bending vibration of a multi-span beam in the x-direction that has been discussed in the 
previous paper[9]. 

It should be noted that if the wave components f^ y ) are replaced by 


</,) - . 

= J“e,x + 2 h n 

Mii = f*e,y + 2 Jy K » ( jy = ) 


( 25 ) 


the dispersion equation, Eq. (23), remains unchanged. For free vibration analysis, we 
may, without loss of generality, choose fi Bjc and as follows 


-» < ^e,x * * > 
~ n < Pe,y * 31 • 


( 26 ) 


The range of and ju^ y specified in Eq. (26) is known as the first Brillouin zone 
(Brillouin 1953) in the case of two dimensional periodic structure. A positive 
corresponds to a wave propagation in the positive x-direction, and a negative fj^ x 
corresponds to that in the negative x-direction. The meaning of is analogous, except 
that it pertains to the y direction. 

The natural frequencies of the grillage can be readily determined from the 
dispersion equation (23), given the values of wave constants (jiq x , /t 0y ), where and r) e 
are both functions of to as seen from Eq. (15). Analogous to the one dimensional case, 
it can be seen from the dispersion equation that the natural frequency of the grillage is 
a multi-valued function of the wave constant. The relationship between X e (or r| e ) and 
Me,y) i n th e propagating bands may be represented by a set of frequency surfaces as 
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shown in Fig. 2. Among four edges of each surface, the edge along which the value of 
Xq (or r| e ) is a constant as well as a maximum corresponds to an upper-bound frequency. 
These upper-bound frequencies can also be determined from one of the conditions, each 
of which is obtained by letting a denominator in the dispersion equation be zero, namely, 


1 - coshX. e cos\ 0 = 0 
simq e = 0 


( 27 ) 


The upper-bound frequencies along the edges - 0 and n, determined from the first 
condition in (27), correspond to a single span beam with clamped ends for bending 
motion. The upper-bound frequencies along the edges of p ey - 0 are n, given by the 
second condition of Eq. (27), correspond to a single span beam with clamped ends for 
torsional motion. It is interesting to note that the propagation band (or frequency surface) 
between 1.5it and 2n for the grillage structure is a non-propagating band for a one- 
dimensional periodic multi-span beam, as shown in [9]. The range of a non-propagating 
band for the grillage structure becomes narrower due to the contribution of the torsional 
motion in the beams. 

Analogously, the dispersion equation for the wave motion of ^(f) has the same 
form as Eq. (23). It may be obtained easily from Eq. (23) by formally replacing the 
subscripts x and © by y and l P, respectively, to yield 
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( 28 ) 


2 Y y 


sinh\, ( cos\, - cosfi^ ) - sinX v (coshX v - cos y ) 


COSX^ coshX^ - 1 


2 / 1 -Yy n v 


COST],,, - COS \/Xy 

sinri,,, 


+ v = 0 
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where 





{ m L* ' 

1 

7 

ml \ 

— 

y y a> 2 


x a> 2 


D 

> Itjr 

C 


r ^ 


* / 


m <yK 


(29) 


v — k 

y y 
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\ 

Dl 

c 2 

X 

7 D 

V 

D 2 

y 

cl) 

L y 2 Ll 

\ y * ) 

k 

ii 

k' 

Ll Ll 
y x ) 


l 

'7 


0 sv s 1 

'y 


(30) 


23 Deflections of Bending and Torsional Motions 

The bending and torsional deflections of the beams in the x and the y directions 
within an individual span can readily be obtained by substituting Eqs. (19) and (20) into 
Eqs. (16) and (17), where © ap is given by Eq. (21). It is noted from Eqs. (19), (20) and 
(22) that the ratios A a ^JB a ^ and C a( JD a ^ are independent of the span numbers a and |3. 
This implies that deflections in two neighboring spans differ only by a phase shift. Thus, 
the deflections of bending and torsional motions of the entire stmcture can be represented 
as follows: 

X a ^) =Z 10 (i)e-'^ (a - 1)+ ^ P1 , ( a, p = 0,±l,±2,...±oc ) (31) 
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where 


«.,(?) - , ( o, p - o,±i,±2,...±» ) 


(32) 


-US) - \ „/>.(?) * a,. 4,(1 -I) < 33 > 

- c », sin %£ * D 0| sinr| o (l-i;) (34) 

and where f k (-) is defined in Eq. (18), and coefficients A l0 , B l0 , C 01 and D 01 are given for 
several cases as follows: 

Case a: fi Q * *j x ii and fi Qj * j y x 

In this case, the wave motions of ©^(f) propagate in four directions, as shown in 
Fig. 3. The coefficients A 10 , B l0 and C 01 and D 0l are obtained from Eqs. (19) to (22) 


Ao -C[<(l)e-^ - /4(0) ] / 

*,o =C[/4(0)e-^ -/4( 1)] 


(35) 


where C is determined 
support (a,P), that is 


C = e r n = 1 

^oi e ’ u o 1 1 


(36) 


by using the continuity condition of angular displacement at the 


i.'Xct) l,„ - ejt) | t „ = 


(37) 


Eq. (37) is solved to yield 


r = L , sinT le 

HK) 


where A(-) is defined in Eq. (20). 


( 38 ) 
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Case b: =y> and fi e<y *j y ii 

It should be noted that when = jji the harmonic wave in the x direction has 
the appearance of a standing wave. As discussed in [9], the deflection in each span in 
the x direction coincides with a mode shape of single span beam with either idealized 
supported or fully clamped ends. The ratio between two angular displacements at two 
ends of a span in the x direction is 

_ Mqp = (-1)'* , s B - integer of [— ] (39) 

©_ ... * 


It follows from substituting Eq. (39) into Eq. (19) that 

afJ 

Thus, the coefficients in Eqs. (33) and (34) are 

*10 - 1 . *10 =(- l ) ,,+1 

C 01 =Ce— , D 0l = C 
where C may be determined from Eq. (37) to yield 

r, ( -iy ^ (1 )‘(- 1)Vl < (0) 

I r sim 1e 


( 40 ) 


( 41 ) 

( 42 ) 


( 43 ) 


Case c: /i Qrt * j and fi ey =j y ji 

Contrary to case b, the wave motion in the y direction has a standing wave 
appearance. The torsional motion in the y direction in each span has the same form as 
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that of a single span beam with either simply supported or fully clamped ends. The 
corresponding angular displacements at the two ends of a span in the y direction are 
related as 

0 n 

■ - — = ( _ l) ,r * 1 , s T = integer part of [_ ®] (44) 

0 a(>-l 71 

which is based on the fact that the odd-numbered (even-numbered) mode shapes of a 
single span beam for torsional motions are symmetric (anti-symmetric) with respect to the 
mid-point of the span for either simply supported or fully clamped ends. Hence, the 
coefficients in Eqs. (33) and (34) are obtained as 

-/4(0)] , 

(45) 

*io =C[/4(0)e-‘>~ -/4(1)] 


in which 


C n , = 1 


01 


An =(-l r 




(46) 


C 


(-iy> 


A sint le 

*(K) 


(47) 


Case d: /x e ^ and fi e<y =jji 

In this specific case, the motion has the form of a standing wave in both the x and 
y directions. Combining the cases b and c, we have 


Aio = L x sinr\ e (-iy- , 


Ao =^ 10 (-1 ) f * +l 


(48) 
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c 01 - [<(1) - (-1)"*‘<(0) ](-!/’ , Z > 01 = C„,(-lf‘ 


(49) 


3. Free Vibration of a Finite-Size Grillage 

3.1 Boundary Conditions along the Exterior Edges 

For a grillage of finite size, having N x spans in the x direction and N y spans in the 
y direction, the boundary conditions along the four exterior edges should be imposed. 
These four edges are identified as a = 0, a = #,, (J = 0 and (J = N y . Each boundary is 
treated as an edge-beam, with one-half the elastic and inertia properties of an interior 
beam. The boundary conditions at each node along each edge may be expressed in the 
following general form 

i - is,, - is„ , ) [( i - - < i - #„.)-*."*( i,o] 

2 1 > ( 50 ) 

where b - is the Kronecker delta. The case k afi = oo represents a clamped edge, for which 
the boundary condition (50) reduces to © a p(0 = 0 (or 0 a p =0). The case 

/ 

^ 4i r-O.N, I-0.AT, 

will be referred to as an ideally supported edge. 

3.2 Mode Shapes of Bending and Torsional Motions 

For a grillage of finite size there generally exist four wave motions in different 
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directions, as shown in Fig. 3, sharing the same frequency to, due to reflections at the 
boundaries. By taking into account these four wave motions, the angular displacement 
function at node (a,p) becomes 


© ap = + c 2 e vv,P) 

(52) 

+ c 3 e + c c *Ov.°**»,P) 

where c t {£=1,2,3, 4} are constants to be determined by imposing the boundary conditions. 
In view of Eqs. (16) through (20), the mode shapes of such a grillage can be obtained by 
combining the above four wave motions. Analogous to the case of a grillage of infinite 
size, these mode shapes are formulated for the following four cases: 

Case a: ^ * jji and * jji 

In this case, four wave motions in (50) propagate in directions not parallel to the 
* or the y axis shown in Fig. 3. Once the value of 0 af> is determined, the mode shapes 
for the bending and torsional motions are obtained from Eqs. (16) and (17) with the 
coefficients given in Eqs. (19) and (20), respectively. 

The case of a grillage with clamped edges along the four edges is selected as an 
illustrative example. The other cases are listed in Table 1. Upon imposing the boundary 
conditions along the edges of a = 0 and 0 = 0, i.e. 


©op = 0 > ( P = 0,1,2,...,# ) 

© a0 =0, ( a = 0,1,2,...,#) 


we obtain C 1 -C 2 --C 3 --C 4 = -C. Thus, Eq. (52) reduces to 


( 53 ) 
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( 54 ) 


©a,p = 4C sin )“e,x a sil V*e,,P 
By imposing on equation (54) the other two boundary conditions along the edges a = N x 
and (5 = N y 


©,.,=0, (P =0,1,2 tf) 

e oAI =0, ( a =0,1,2„..,JV, ) 

we then have 

sil V*e,A “ 0 ; sir Ve,y N y = 0 


Therefore, 


(55) 


(56) 


^e,x 



(; - 1 , 2 , 


A-i) 


(57) 


^ (* = l,2,...,tf -1) 

N 

y 

Table 1 lists the angular displacement functions @ ap and the equations for obtaining the 
wave constants and p e for all sixteen sets of boundary conditions. 

It is seen from the Table 1 that when a pair of opposite edges of a grillage are 
either ideally supported or clamped, the associated wave constants in that direction can 
also be an integer multiple of Jt. These cases are discussed in detail in the following 
sections. 

Case b : ^ =jpi or fi e<y *jyn 

In this case, a pair of opposite boundary edges parallel to the y direction can be 
either ideally supported or clamped. Then the @ a p(f) motion can be represented by two 
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waves propagating in the positive and negative ^-direction. In view of Eqs. (31) and (32), 
the mode shapes for the bending and torsional motions may be written as follows 

*„,(?) - * V <"•>') (58) 

- {c.o. 1 (5)«'* wa, '* ) - c J e; i (a«r^<'-‘»}(-iy-" (5») 

where asterisk denotes the complex conjugate and 0 ( § ) and 9 0l (^) are given in Eqs. 

(33) and (34) with the coefficients specified by Eqs. (41) through (43). The coefficients 
c i an d c 2 * n Eqs- (58) and (59) are determined by using the boundary condition along 

edge P = 0 for the following two cases, namely, the ideally supported and the clamped 
edges. 

Consider first the case in which the boundary along 0 = 0 is clamped, i.e. 
®„o = ^aoC 1 ) = Th e opposite boundary along p = N y could be either the ideally 
supported or clamped. Imposing the boundary condition © a0 = 0 on Eq. (59) results in 
c i = - c 2 . Thus, the mode shapes for the bending and torsional motions reduce to 


*,(1) *Jr l0 (?)(-iy- ( “-‘>si^p 

(60) 


(61) 


where Im{-} denotes the imaginary part of a complex expression, and the common factor 
2 ic 2 has been omitted. 

Consider next the case of ideally supported edge along p = 0. The opposite 
boundary along P = N y may be either the ideally supported or the clamped edges. The 
boundary condition along p = 0 obtained from (50) has the form of 
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( 62 ) 


4*.®.. - J0,[-C„(O) -*".(!)] * C 61,(0) , 


( a = l,2,-,N x -l ) 


Substituting Eqs. (58) and (59) into Eq. (62), we conclude that if c t = c 2 . Then Eq.(62) 
is equivalent to the dispersion equation Eq. (23) in the special case of = jpi. It 
follows from Eqs. (58) and (59) that the mode shapes for the bending and torsional 
motions are 



( 63 ) 


( 64 ) 


where Re{-} denotes the real part of a complex expression. 

Of course, Eqs. (63) and (64) are equivalent to Eqs. (60) and (61) when one edge 
is ideally supported and the opposite edge is clamped. 

Case c: *jpi and =j y i t 

Analogous to Case b, the motion can be represented by two waves propagating in 
the positive and the negative .^-direction. Two opposite edges parallel to the x direction 
can be either clamped or ideally supported. The following mode shapes of the bending 
and torsional motions are obtained from Eqs. (31) through (34), with the coefficients 
given by Eqs. (45), (46) and (47) 


-Ml) - c.-W?)* 


•‘A./ 0-1 ) 


C 2 ^1 0' 






(-ly '’ 1 


( 65 ) 
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K,K) - e cl (t) ( -1 /-""’[c.e * c 2 e^-‘ } «*> 

Let the edge along a = 0 be clamped; i.e. @ op =0, where p = l,2,...fl y -l. The opposite 
edge along a = N x may be either clamped or ideally supported. By imposing the 
boundary condition along a = 0 on Eq. (66), we find c x = - c 2 . Thus, the expressions for 
mode shapes Eqs. (65) and (66) reduce to - 



( 67 ) 

9„,(« =e 01 (S)(-iy.« | -‘W 9; ,a 

( 68 ) 


If the edge along a = 0 is ideally supported, then the boundary condition may be 
obtained from Eq. (50) as follows 


\K ©o, = ^X(O) + 0o P -t (O) - ®o P (l) ] 


( 69 ) 


( p = 1,2, ...,1V -1 ) 


The edge along a = N x may be either ideally supported or clamped. Substitution of Eqs. 
(65) and (66) into the boundary condition Eq. (69), and use of dispersion equation (23) 
result in c l = c 2 . Hence, the mode shapes in this case become 


*.,<?) = «e{Ar i0 (l) e -l/.' 

5„,(« =8„,(«(-iy- < '“We..a 


( 70 ) 

(71) 


Case d: =y> and = j y n 

In this case, each pair of opposite edges must be both clamped or both ideally 
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supported, and the mode shapes coincide with the deflection shapes of a grillage of 
infinite size. These mode shapes may be expressed as follows: 



(72) 


(73) 


where X 10 (|) and 0 O1 (£) are given by Eqs. (33) and (34), and the associated coefficients 
are specified by Eqs. (48) and (49). 

33 Example and Discussion 

As an example, a grillage of finite size, consisting of four spans each in the x and 
y directions (N x = 4 and N y = 4) with all clamped edges is investigated. The 
nondimensional parameters are chosen as v, = 0, y* = 0.5, k* = 0.08. 

The wave constants and the associated frequencies obtained from the first and 
second frequency surfaces are listed in Table 2. It is seen from Figs. 2(a) and 2(b) that 
the first upper-bound frequency occurs at = 0, whereas the second upper-bound 
frequency appears at = n. Note that the first upper-bound frequency corresponds to 
a bending motion of the beams in the x- direction, whereas the second upper-bound 
frequency corresponds to a torsional motion of the beams in the y-direction. For 
emphasis, each of these is marked by an asterisk in Table 2. 

The mode shapes in the spans in the x- and y-directions are shown in Figs. 3(a)- 
3(m) for both bending and torsional motions, where the parameters are fixed at v* = 0.5, 
Y* = 0.5, k,. = 0.08. It is seen that the mode shapes of the bending and torsional motions 
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for all spans are related, since they must be consistent at the nodes. Moreover, the mode 
shapes must be either symmetric or anti-symmetric since the boundary conditions on 
opposite edges are symmetric in this case. The mode shape corresponding to each upper- 
bound frequency in the first column in Table 2 is devoid of torsional motion. In contrast, 
the mode shapes corresponding to each upper-bound frequency in the second column is 
devoid of bending motion. 

4. Conclusions 

Exact solutions have been obtained for wave propagation in a periodic grillage of 
infinite size, and for free vibration of one of finite size. Attention is focused on the 
motion at the periodic supports, thus by-passing the complicated motion between two 
neighboring supports. It is shown that the non-propagation frequency bands are narrower 
for a two-dimensional grillage structure than that of a one-dimensional periodic multi-span 
beam. The mode shapes of free vibration for both bending and torsional motions obtained 
herein may be used in a forced vibration analysis using the normal mode approach. One 
implicit assumption in the present analysis is that wave propagation in the x and y 
directions are synchronized such that they appear as components of a "plane" wave 
propagating in an oblique direction. The possibility of a "curvilinear" wave front will not 


be considered at this time. 



Table 1. Angular displacement and wave constants of grillage 
of finite size with various boundary conditions. 


Boundary conditions 

MQ,y* m yn) 

Equations for and fx^ y 

Cx~Cx an< ^ Cy'Cy 

sin^^a siiyr^P 

sin/ie^V x =0 ; sin 1 u e yA/y=0 

S x -S x and C y -C y 

cos/^a sin^p 

sin^e^V^O ; sin^/v^O 

S x -C x and C y -C y 

cos^a sin^e y P 

cos/i 0 ^Af x =O ; sin - u ey A/y =0 

C x -S x and C y -C y 

siiyr^a siry^ 

cos / u e ^V x =0 ; sin^eyNy=0 

C X C X and S y -S y 

siiy^a cos/^ y p 

sin/i e ^V x =0 ; sin^QyA/y =0 

C x -S x and S y -S y 

siry^a cos/ieyP 

cos/^jN^O ; sin^//y =0 

S x ~Cx an< i S y -S y 

cos^a cos/^yP 

cos / u e ^A r x =0 ; sin l a e yAry =0 

S x -S x and S y -S y 

cos/^ct cos/^yP 

; sir^ ey N=0 

C x -C x and C y -S y 

sin^e^a sin/r^p 

sin/JvJf=0 ; cos/^yA/y =0 

C x C x an ^ 

sin^e^ct cos^p 

sin^^V x =0 ; cos i u Qy A r y =0 

C x -S x and C y -S y 

siiy^ct sin/i^P 

cos l « e ^A r t =0 ; cos / u e yATy=0 

C x -S x and S y -C y 

sin^e^a cos^p 

cosf^Jf=0 ; cos J u e yA ^=0 

S x -S x and C y -S y 

coSyUe^a siry^yP 

sin^V =0 ; cos ^/^0 

S x -S x and S y -C y 

cos^e^a cos/^yP 

sinyUe^Ar x =0 ; cos^yA^O 

S x -C x and C y -S y 

cos^a sin^eyp 

cosyu e ^ r =0 ; cos^eyA/y =0 

S x -C x and S y -C y 

cos^a cos^yp 

cos/ie^^O ; cos j a e> /^=0 


Note: C x> C y = clamped support; S x , S y = ideal support; subscripts indicate direction. 
M x = number of spans in the x-direction; N y = number of spans in the y-direction. 
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Table 2. The first two groups of frequencies of a finite-size grillage( N x = 4, N y =4 ) 
with all edges clamped ( v x =0, y x =0.5, q x =0.08 ) 


in in 

** m TT 

* y 

* J 

c 

Jt 

* 4 2) 

n 

0 1 

1.5056* 

— 

0 2 

1.5056* 

— 

0 3 

1.5056* 

— 

0 4 

— 

— 

1 1 

1.3829 

1.6852 

1 2 

1.4198 

1.7664 

1 3 

1.4422 

1.8882 

1 4 

— 

1.9947* (ri e =Jt) 

2 1 

1.2162 

1.7729 

2 2 

1.2752 

1.8386 

2 3 

1.3185 

1.9332 

2 4 

— 

1.9947* (rj e =n) 

3 1 

1.0666 

1.8207 

3 2 

1.1478 

1.8775 

3 3 

1.2076 

1.9531 

3 4 

— 

1.9947’ (r| e =n) 


Note: An asterisk denotes an upper-bound frequency. 
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Fig. 1(a) Periodic grillage with rotational constraints at each node 



Fig. 1(b) Notations for displacements at node (ct,P) 




Jl 


(a). The first frequency surface 



(b). The second frequency surface 
Fig. 2 Dispersion relationship of wave motion in a grillage 
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Fig. 3 Wave propagation in four directions in a grillage of finite size 
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bending mode 


torsional mode 



4(b) Me* = 3ji/4, = 2ji/4; Xq = 1.1696ji 
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bending mode 


torsional mode 



4(d) = 2it/4, jx Qy = ji/4; Xq = 1.2369ir 
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4(f). = 2 ji/ 4, //Qy = 3it/4; Xq = 1.3288jt 
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4(g)- fas - n/A > fa.y = n/4; ^9 = 1.394631 


bending mode 



4 ( h )- fas = x/4, fiQ y = 2jt/4: 


torsional mode 



= 1.426631 
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4(0- Mu* = «/4, Mb, y = 3jc/4; Xq = 1.4462ji 



4(j). ^ = 0, //e.y = 3c/4; Xq = 1.5056jt 
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4(k). fa^ - 0, fa y = 2^/4; = 1.5056ji 



4(1). fa^ = 0, fa y = 3n/4; = 1.5056n 

Fig. 4 Mode shapes of grillage with clamped edges ( \ x - y x = 0.5; k x = 0.08 ) 
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Abstract 


The paper deals with random vibrations of the space shuttle weather protection systems. 
It is shown that the Timoshenko beam theory must be utilized to describe the structural behavior 
of the system. Use of the simple Bemoulli-Euler theory may result in an error of about 50% in 
determining the mean-square value of the bending moment in the weather protection system. 


Introduction 


This study deals with the random vibrations of space shuttle weather protection systems 
modeled as Timoshenko beams. The experimentally determined excitation, rather than assumed 
expression of its cross-spectral density will be utilized. Random vibrations of Timoshenko beams 
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under simplest, "rain-under-roof" excitation, namely with both time-wise and space-wise white 
noise, was considered previously. Samuels and Eringen [1] were the first authors to study 
random vibrations of Timoshenko beams. They concluded, for the above excitation, that the 
results produced by using either the Timoshenko or Bemoulli-Euler theory differed by less than 
five percent. Crandall and Yildiz [2] studied effects of both the dynamical models utilized and 
of the postulated damping mechanism. They considered, as an excitation, the time-wise band- 
limited white noise and investigated the growth pattern of the response characteristics with the 
increase of the cut-off frequency. Banerjee and Kennedy (1985) investigated the effect of axial 
compression, whereas Singh and Abdelnaser (1993), studied the effect of boundary conditions. 
The latter study is in agreement with the paper by Samuels and Eringen (1957) in the sense that 
for the length-to-depth ration Lid above nine, the percentage-wise difference between results 
produced on one hand, by the Timoshenko theory and, on the other hand by the Bemoulli-Euler 
theory is less than five percent. Therefore, one may conclude that for this case of the length-to- 
depth ratio one should utilize a simpler Bemoulli-Euler theory. 

However the above studies were concentrating on the excitation whose spectral density 
is constant. As a result, within the modal analysis approach, only the low end of the frequency 
spectrum is contributing significantly to the response. However, the shear deformation and rotary 
inertia are affecting considerably only the higher frequencies. Since the contribution of jth mode 
is inverse proportional to to * for the viscously-damped beam, the significantly affected 
frequencies contributed very little in formulation of the response. This is the reason of 
uncovering a small numerical difference between the Bemoulli-Euler and Timoshenko beams. 
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However, the situation must not be such in the general case of excitation, as was first 
demonstrated by Elishakoff and Lubliner (1985). They considered the band-limited white nose 
excitation with two out-off frequencies, namely, the lower one (n cl and the upper one o) c2 , i.e. 
the spectral density was taken constant for frequency ranges - s w s - (o < . 1 

and to s (o s o) Such a band limited white noise tends to ideal white noise when 

tends to zero and tends to infinity. It was demonstrated that when oj^ is zero or close to 
zero, both the Timoshenko theory and the Bemoulli-Euler theory produce similar results. 
However, when co <;l increases, the higher modes become more important. Since the higher 
natural frequencies are reduced by the effects of shear deformation and rotary inertia, the 
response predicted by Timoshenko theory is in excess of that predicted by the Bemoulli-Euler 
theory. In such circumstances one must utilize Timoshenko beam theory. In several example 
cases Lubliner and Elishakoff (1985) demonstrated that the application of the Bemoulli-Euler 
theory would yield up error of order of 50% or more. 

In this study we consider the response of weather protection systems to random 
excitations. The set of consistent differential equation, as discussed by Elishakoff and Lubliner 
(1985) is utilized. Namely, the fourth order derivative of the displacement with respect to time, 
appearing in the original Timoshenko equation is neglected as discussed by Love (1927), Tseitlin 
(1961), Egle (1969), Elishakoff and Abramovich (1992). 

The typical weather protection systems for launch vehicles are made of thin corrugated 
metal sheets (see Fig. 1), with a much larger bending stiffness in the direction perpendicular to 
the corrugation, namely, I z >>/ y ;. These sheets are supported by thin beams parallel to the O z 
direction. Such a weather protection system is susceptible to excessive shear deformation. It is 
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therefore reasonable to model such a system as a Timoshenko beam. 


Equivalent Timoshenko Beam 

For the sake of analysis, a typical segment of the corrugated sheet may be substituted 
by an equivalent I beam as shown in (Fig. 2). We equate an infinitesimal element of the web 
of the corrugated beam with its equivalent part on the / beam, (Fig.2). 


y - T]cosa 


t { dy = tdr\ ; 


t 

cos a 


( 1 ) 

( 2 ) 


Furthermore, the shear flow in the two elements must be equal, 

xt = x l t l (3) 

where x and x t are the shear stresses in the corrugated sheet and the I beams, respectively. 
Hence 


The shear strains are given by 


Xj = xcosa 


(4) 


Y 




(5) 


where G is the shear modulus of the material in the corrugated beam, and G v is the equivalent 
shear modulus in the I beam. 

Now, for the I beam web to have the same shear stiffness, the displacement u in the 
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longitudinal direction must be the same as that of the corrugated web, 


du = ydr\ = y v dy 


Hence 


and 



cosa 


T i _ xcosa _ x 
G 1 Gj Geos a 

From Eq.(8), one obtains the equivalent shear modulus 


G v = Geos 2 a 

one can show that shear stress energies are equal in the two models 


dU 


1 

2 G[ 


t^dy 


1 x 2 cos 2 a 

2 Geos 2 a 


tdr\ 


1 x 2 


2 G 


tdr\ 


( 6 ) 


( 7 ) 


( 8 ) 


( 9 ) 


( 10 ) 


Free Vibration of Timoshenko Beam 

The equation of motion for a Timoshenko beam may be written as follows (Elishakoff 
and Lubliner 1985), 
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J-.J d*w . d 2 w T 

El + pA - ol 

dx A dt 2 


1 +- 


kG, 


3*w 


dx 2 dt 2 




. q dw p IE d 3 w , \ 

+ P A fW ~ . TT = ^>0 


/ 


-Elfl + pl^l 
dx 2 dt 2 


\ 


( 11 ) 


0 dt kG x dx 2 dt J ' 
where (3 0 is the transverse viscous damping coefficient p o = c/pA. The natural frequencies of 


the structure are found by letting p o ■ 0, <l(x,t) m 0 and 

Mx,t) - ■^(x)e w>t 

For a simply supported beam, the mode shape can be chosen as 

imc 


( 12 ) 


(13) 


%(*) = sin - L 

where j is the number of half sine waves in the x direction. The natural frequencies are 
computed from 


co, = 


El ( jn \ 4 


(j*) 2 

t \r 

i. ^ 

pA( L j 

A 

UJ 

k ' G >)\ 


(14) 


The shear correction factor k is given by (Cowper, 1966) 

k = 10(1 +v) (1+3 m) 2 /B 
where v is Poisson’s ratio, and B is defined as 


(15) 


B = 12 + 72m + 150m 2 + 90m 3 + v(ll + 66m + 135m 2 + 90m 3 ) 
+ 30/1 2 (m + m 2 ) + 5v/i 2 (8m + 9m 2 ) 


(16) 
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Random Vibration of Timoshenko Beam 

Consider a random acoustic excitation field characterized by a cross-spectral density 

<& q (x ir r 2 ;u>) = S g (u))exp[- a(a))|x, - x 2 \ - ry( a>)(x 1 - x^] (18) 

where S q (u>) is the spectral density at the reference point x 2 = x 2 = x; a(to) and y( w ) are > 
respectively, the decay and phase functions of co. They are determined from experimental data 
and are given in the Appendix . Assuming that the beam has reached the state of probabilistic 
stationarity, the cross-spectral densities of response for displacement and moment are respectively 
found as follows (Lin 1976; Elishakoff 1983): 


0 DD (x ir x 2 ;co) = b 2 L ^9 f (y v y 2 ; u)H D (x v y { ,u>)HZ{x v y„ u)d yi dy 2 


(19) 


®Mkk X V X V W ) = b 2 Jo 


( 20 ) 


where H D and H M are the frequency response functions of displacement and bending moment, 
respectively, L is the length of the beam, b is the width of the beam. 

In order to derive the function H D , let the harmonic point loading be imposed at location 
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q{x,t) = 6(* - y)e^ 

Thus, the responses of displacement and moment at location of x will be 

Mx,t) = H D (x,y;oi)e “* 


( 21 ) 


( 22 ) 


M(x, t) = H M (x,y;o})e‘ ait 


(23) 


Substituting Eq.(22) into Eq.(ll) and expending H D in the series in terms of the mode shapes 
\ pjx), we obtain the following expression for H D (x,y;m) 

oo 

HJx,yM • Y. (24) 

with Vj 2 being the norm of jth mode namely, 

v? = J^y (x)dx = £ ( 25 ) 

The modal frequency response function in the j- th mode Hj(to), is derived as 

(26) 

with two nondimensional parameters 


Hf co) = 


1 +aX 1 X 2 [(jJi) 2 - p L 2 /£a> 2 ] 


pA[l+aX,(l + \A(jn)‘ 


a > 2 - to 2 + iu>fJ 0 


1 +aX 1 X 2 (/Ji) 2 
l+aX^l+X^O'n) 2 


X 


i 



(27) 


Moreover, a is an artificial parameter, with a = 1 corresponding to the Bresse-Timoshenko 
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theory and a = 0 associated with the Bemoulli-Euler theory. 

Note that the relationship between the bending moment M(x,t) and displacement w(x,t) 
is 


M(x,t) = -El 


w\x,t) - -2—V ( x,t ) 
V ' kGA y ' 


(28) 


= - EI\ 


W'\x,t) - _J_[c 0 vv(x,t) + pAMx,t) - q(y,t ) 

kGA 


Substituting Eqs.(21)-(23) into equation Eq.(28), we obtain H u as follows 

co 2 - wQJHJwj*) + ^( x " >oJ < 29 > 

The cross-spectral densities for displacement and moment can be rewritten by substituting Eqs. 
(24) (29) into Eqs. (19), (20) as follows 


®dd( x v x 2’^) = b 2 £ E 

y« i 


(30) 
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00 00 

i***®) = b2 ( Er > 2 u2 £ [^"(^^*"(^ 2 ) + nV + 

;*l **1 

+ri(a ) 2 + io)P 0 )Ti» ; .''(x 1 )^ i (x 2 ) + rj(co 2 - 

+ -^-£ [n(t*> 2 + iojPo)^*^ + %\xJ]Hfo)J?(x l ;ia) 

pA k 

+ Jl-53 [ T l( a>2 “ i^Po)^W + V'fX^fvVfe™) 

pA j 


( 31 ) 


JjL* 

pA 


fc/W") 


ap 

ri = —L. 
kG 


where two integrals of I jk and J t are defined as follows: 


//to) = 


v,v t 


(32) 


and 


Jj(x) - £ f L &{y*,viNfy)<ty 

\ 1 . j® 


Note that the integral of //to) possesses the Hermitian property. 


(33) 


//to) = 7/<o) 
// _ to) = /;(o>) 


(34) 


Hence, the mean square values of both displacement and moment can be obtained by integration 
of both two spectral densities defined by Eqs (28), (29) over positive range of frequencies only, 
i.e. 
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(35) 


E[w\x,t )] = 2^ <& DD {xjc,u>)dw> 


E[M\x,t )] = 2jJ<Z> MM (xj:,i>i)dw 


(36) 


in which, the integrals l )k and J j can be replaced by their respective real parts and are evaluated 
in the exact form as follows: 


l ^) = -xr i "l L Re {®ftvyv™)}y{y^{yd d yv d y2 

1 + exp[-(a + ry)] 


V;V* 


* 45 o (co);7ot 2 


\2Re\ 


[(a + ff) 2 + (foi) 2 ][a + ry) 2 + (jn ) 2 ] 


+ab 


l+[a 2 + pJO'ibi 2 )' 1 


r rz 


[(a + ry) 2 + (jn) 2 ][(a - ry) 2 + (far) 2 ] 


(37) 


7(*,u>) = _i.J^ L /?e{O g (y^;(o)}^(y)tfy 


; 


= 2S (oi)Re L. _1 ——[(a - ry)sin(/Ji£) - JJicos (/ji£) + ~ (a " ,7)l ] 

[(a - ty) 2 + (jn) 2 

+ __ i [(a + iy)sin(/jr^) + jncos(jn%) + (-1)/Jie (a + ^ * ' l >]| , 

(a + iy) 2 + (jn) 2 j 


(38) 


where nondimensional parameters a and y are defined as follows: 


a = a(u>)£ ; y = y (a>)L ; | = x!L 


(39) 
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Numerical Example and Discussion 


A thin corrugated metal sheet-the typical weather protection system for launch vehicles, 
(Fig.l) is chosen as an example for application. We chose an element of a sheet ( Fig. 2) and 
modeled it as an equivalent /-beam, with properties / = 0.0247 in 4 , A = 0.264 in 2 , b = 1.1 in. 
Thus, two nondimensional parameters in Eq.(27), become Xj = 8.852; Xj = 3.3 x 10 s . 

Table 1 lists the natural frequencies of the beam within both the Timoshenko theory, and 
the Bemoulli-Euler theory approximations. It can be seen, as expected, that the values of natural 
frequencies of the Timoshenko beam are less than those of the Bernoulli beam , and for the 
higher order modes there is a bigger difference between two approximations. 

Fig. 4 and Fig. 5 portray the exact spectral densities of the displacement and the moment 
at the mid-section of the beam. The solid curve presents results associated with the Bemoulli- 
Euler theory while the dashed one denotes the results obtained using the Timoshenko beam 
theory. Only odd-numbered modes contribute to the response, due to sinusoidal mode shapes. 
It is shown that within the same range of frequencies (0 to 1000 Hz), the spectral density of 
response for the Timoshenko beam has more peaks than that for the Bemoulli-Euler beam. In 
other words, for the Timoshenko beam there are generally more modal contributions to the 
response of the structure due to the increased modal density of natural frequencies. For the 
spectral density of displacement ( Fig. 4) the first peak is the largest within both beam theories. 
The difference of values of the first peak for two beams is very small. Hence the value of areas 
under the spectral density curves, or the mean square values of displacement are very close, with 
attendant difference constituting about 1.2%. 

However, for the spectral density of moment (Fig. 5), the higher order modes dominate 
in formulating the response of the structure. This is because there is a high peak around 585 
Hz in the loading spectral density S q ( to) (Fig. 3), and in addition there appears a factor fk 2 in 


5-12 




the expression for H M . Therefore the contribution of higher modes is significant. The 
difference between the values of spectral densities calculated via the Timoshenko and the 
Bemoulli-Euler theories, and their associated contribution to the mean square value of moment 
is remarkable. 

The value of mean square for the Timoshenko beam is almost twice that for the Bernoulli 
beam. This implies that choosing Bemoulli-Euler beam model may not be on the safe side, for 
the design of the weather protection systems. As a result the refined Timoshenko beam model 
rather than simple Bemoulli-Euler theory should be chosen for the space shuttle weather 
protection structure under the acoustic loading. 
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Table 1. Natural Frequencies of the Weather Protection System Modeled as 
Bernoulli Euler-Beam or Timoshenko-Beam 



Natural Frequencies, f ; [Hz] 

i 


Timoshenko-Beam 

Bernoulli Euler-Beam 

i 

1.739829 

1.742619 

2 

6.926167 

6.970474 

3 

15.46189 

15.68357 

4 

27.19260 

27.88190 

5 

41.91663 

43.56546 

6 

5 9.39734 

62.73426 

7 

79.37590 

85.38831 

8 

101.5836 

111.5276 

9 

125.7528 

141.1521 

10 

151.6251 

174.2619 

11 

178.9584 

210.8569 

12 

207.5304 

250.9370 

13 

237.1413 

294.5025 

14 

267.6140 

341.5533 

15 

298.7939 

392.0891 

16 

330.5479 

446.1104 

17 

362.7619 

503.6167 

18 

395.3398 

564.6084 

19 

428.2004 

629.0853 

20 

461.2764 

697.0474 

21 

494.5116 

768.4949 

22 

527.8599 

843.4274 

23 

561.2834 

921.8451 

24 

594.7518 

1003.748 

25 

628.2399 

1089.137 

26 

661.7276 

1178.010 

27 

695.1992 

1270.369 

28 

728.6419 

1366.213 

29 

762.0460 

1465.542 

30 

795.4041 

1568.356 

31 

828.7104 

1674.656 

32 

861.9608 

1784.441 

33 

895.1523 

1897.712 

34 

929.2831 

2014.467 

35 

961.3519 

2134.708 

36 

994.3586 

2258.434 

37 

1027.303 

2385.645 

38 

1060.186 

2516.341 

39 

1093.008 

2650.523 

40 

1125.770 

2788.190 
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Appendix 


The model of cross-spectral density for acoustic loading utilized in this study is given 
by Eq. (18) where the spectral density S/co) is approximated by the following expression 


= 


3 


E 


CO 


— 

exp 

CO 


} 



b f - (i - <^ 2 

2(1 ~c) 


u). = 2nfj ; b. = 


CO 

CO. 
\ > 



(40) 


with parameters listed in the Table: 


j 

1 

2 

3 

a i [ ( psi) 2 /(r/s)] 

1.1 x 10' 2 

4.5 x 10' 3 

1.9 x 10’ 2 

f ;■ [Hz] 

105 

290 

585 

S 

0.8 

-1.0 

0.995 


(41) 


(42) 


Decay function a(to) is defined as 

a(u>) = 0.03 + 5.0 x 10' 5 -^-f 
and has a width dimension of [1/ft]. Phase function y(to)reads 

y(a>) = 0.77 sin (0.14-il) + 0.045_^ , 

2ji 2n 

having a dimension [degree/ft]. 
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ij)(Ax,a)) S q (to) [(psi) 2 /(r/s)] 



cl»/2jt: [Hz] 

Fig. 3a Spectral density of acoustic loading 
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MM 



Fig. 4 Spectral densities of the displacement at the mid-point of the beam 



co/2jt 


Fig. 5 Spectral densities of the bending moment at the mid-point of the 


beam 
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ABSTRACT 

Several recent developments in the stochastic linearization technique are summarized in this review paper. The 
nonlinear oscillator subjected to colored noise is examined; the case of the nonlinear damping is discussed; the "true" 
stochastic linearization technique is described. The results of the latter match the exact mean square responses of 
the nonlinear structure. The combination of the stochastic linearization with the Monte Carlo method is outlined. 
In addition, the accuracy of a new linearization technique in contrast with the classical linearization scheme is 
examined for a Duffing oscillator subjected to white or colored noise excitations. The results obtained by the two 
linearization schemes are compared in terms of percentage-wise error in reference with the exact solution or 
numerical results obtained through Monte Carlo simulation. These applications confirm a superior performance of 
new linearization technique in comparison with the classical one, in several examples considered. Under some 
circumstances however, namely for some nonlinear softening oscillators, the conventional linearization may yield 
more accurate results. The method of weighing functions improves the accuracy of both the conventional and new 
stochastic linearization methods. 

The developments which are described in this review mostly took place after excellent accounts on the classical 
version of the stochastic linearization technique, the monograph by Roberts and Spanos (1990), and the review article 
by Socha and Soong (1991), have been published recently. 


INTRODUCTION 

Stochastic linearization technique is the most versatile method for analysis of general non-linear 
systems and structures under random excitation. In almost thirty years since it was proposed for 
the first time by Booton (1953) and Kazakov (1954) it has been widely applied for the study 
of various non-linear systems which are not amenable to exact solutions. For example, the 
monograph by Roberts and Spanos (1990), and the recent review paper by Socha and Soong 
(1991) give comprehensive accounts of some of these developments. The fundamental idea of 
the method lies in replacing the original non-linear system by a linearized one in such a way that 
the difference between two systems is minimal in some probabilistic sense. Following the 
classical approach, the linearized system parameters are determined in the manner that the 
difference between the nonlinear force and the force in the linearized system is minimal in the 
mean square value. In some recent papers by Falsone (1992.a,1992.b) it has been demonstrated 
that, when parametric excitations are present, it is more suitable to measure the difference on the 
coefficients of the Ito differential rule. It is remarkable that, when only purely external 
excitations are present, Falsone’s approach coincides with the classical stochastic linearization. 
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Recently, for the case of purely external excitations, some new stochastic linearization techniques 
have been suggested (Wang and Zhang, 1985; Zhang, 1989; Zhang, Elishakoff, Zhang, 1991; 
Elishakoff, 1991; Elishakoff and Zhang, 1992.a; Elishakoff and Zhang, 1992.b) in which the 
differences between the nonlinear original system and the linearized one are considered in terms 
of potential energy. In particular, Elishakoff (1991) investigated a Duffing oscillator subjected 
to a white noise excitation; it has been shown that the best results are obtained when the 
linearized system parameters are obtained through minimizing the mean square error between the 
potential energies of the original non-linear and the replacing linear system The exact solution 
for the stationary probability density of Duffing oscillator is readily available and approximate 
solutions are not needed for this problem; therefore the validity of the modified stochastic 
linearization technique is easily checked for this case. Moreover, the results furnished by the 
stochastic linearization technique are compared with those yielded by the Monte Carlo 
simulations when the exact solutions are unavailable. 

It has demonstrated, (Elishakoff, 1991) that for some combinations of the parameters of 
the Duffing oscillator the proposed linearization criteria may yield results in perfect agreement 
with exact solution. For other sets of parameters, the proposed linearization yields results which 
are slightly greater than the exact probabilistic responses, whereas the conventional linearization 
yields responses which are below the exact values. Since generally engineers utilize the notion 
of safety factors, the structures designed through use of conventional stochastic linearization, may 
turn out to be "over-designed." When the Duffing oscillator is extremely weakly nonlinear, the 
conventional stochastic linearization may exhibit less error, than the energy-wise stochastic 
linearization. However, for these almost linear systems the percentage-wise error is under two 
percent for either stochastic linearization criteria. Therefore, for the Duffing oscillator, the 
energy-wise stochastic linearization is almost always preferable. 

The aim of this paper is to give a brief account of some recent developments, as well as 
to investigate the accuracy of this new linearization technique in the case of colored noise 
excitation. The results obtained in this way have been compared with those obtained by means 
of the Monte Carlo simulation. 

The combined use of concept of potential energy and energy dissipation function is discussed for 
nonlinearly damped structures with nonlinear restoring force. Finally, the combined use of the 
stochastic linearization and Monte Carlo method is discussed. 


NEW VERSUS CLASSICAL LINEARIZATION TECHNIQUES 

Let us consider the single-degree-of-freedom system governed by the following equation of 
motion 

mX + cX + fiX) = W(r) t 1 ) 


where m is the mass, c is the damping, /(•) is the restoring force, which is considered here as 
a non-linear function of the displacement X, W(t) is a random, stationary, Gaussian excitation 
upper dots denote time derivatives . The basic idea of any stochastic linearization technique 
consists in the replacement of the original non-linear equation (1) by such a linear equation that 
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the difference between the two systems is minimal in some probabilistic sense. According to the 
classical linearization technique (Booton, 1953; Kazakov, 1954; Roberts and Spanos, 1990; Socha 
and Soong, 1991), equation (1) is replaced by the following linear one: 

mX + cX + k®X = W(0 (2) 


in which is the spring coefficient, chosen in such a way that the difference between equation 
(2) and equation (1) attains a minimum in mean square sense, that is 


£p«) - O) 2 


min 


(3) 


where E[-] means mathematical expectation of (■)• Following this approach one finds the well- 
known relationship for k£\ namely (booton, 1953; Kazakov, 1954) 

= Eimx\ (4) 


where E[f(X)X] is evaluated taking into account that the replacing system is linear; so that, since 
the input is Gaussian, the response is assumed to be likewise Gaussian. 

An alternative criterion was suggested by Kazakov (1954) in his pioneering paper. This 
criterion requires that the mean square values of the nonlinear restoring force f(X) and the 
equivalent, linear restoring force k^ 2) X be equal 

£[/ ! tf)] = £[(Crf] (5) 


which leads to expression 



E[f\X)} 


M 


E(X*) 


( 6 ) 


It is interesting to note that the Booton- Kazakov criterion (3) is almost universally utilized in the 
literature. To the best of author’s knowledge, the Kazakov criterion (5) has been elucidated only 
in the book by Popov and Paltov (1960). They have observed that for numerous cases, mean- 
square values of the response, furnished by the Booton- Kazakov criterion constituted a larger 
value than the result delivered by the Kazakov criterion. Hence they suggested to utilize the 
arithmetic mean of these two results to obtain an approximation which would be closer to the 
exact mean square, than the ones delivered by either Booton-Kazakov (1953) or Kazakov (1954) 
criteria. 


Bolotin (1984) demonstrated on the particular example of a half-degree-of- freedom system 
(the simple oscillator with negligible mass, for which an exact solution was constructed by 
Caughey and Dienes, 1961) that the exact mean-square displacement may be in access of results 
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furnished by either Boo ton- Kazakov or Kazakov criteria. Therefore the approximation suggested 
by Popov and Paltov (1960) is not necessarily more exact than those given by either criteria (3) 
or (5). 


Recently, a new stochastic linearization technique has been suggested (Wang and Zhang, 
1985; Zhang, 1989; Zhang, Elishakoff, Zhang, 1991; Elishakoff, 1991; Elishakoff and Zhang, 
1992), based on the concept of potential energy. 

Following this technique, equation (1) is replaced by the following linear one: 

mX + cX + k™X = W(t ) , C 7 ) 


in which the new stiffness spring k^ 2) is chosen in such a way that the mean square deviation 
between the potential energies possessed by the original non-linear system (1) and by its linear 
counterpart in Eq.(5), attains its minimum, that is: 


{u(X) - *® r %} 2 = min 


( 8 ) 


resulting in: 



. 2 E[UW£] 
E[X 4 ] 


( 9 ) 


In equations (8) and (9) the expression U(X) represents the potential energy of the original 
non-linear system. Elishakoff (1991) applied this new technique to a Duffing oscillator subjected 
to a white noise excitation, and performed a systematic comparison with an exact solution. It was 
shown that the new technique exhibits a "better" performance than the classical linearization; 
namely, for particular values of parameters, the new linearization may yield the exact response. 
In the next section the classical and new linearization techniques are applied to Duffing oscillator 
to illustrate the superiority of the new stochastic linearization technique. 

In analogy with Kazakov’s (1954) criterion, given in Eq. (5), Elishakoff and Zhang (1992) 
suggested to use the fourth criterion, based on the requirement that the mean-square values of 
the potential energy of the original system and its linear counterpart be equal: 

eM = (10) 


yielding 



£{[W] 2 ) 


N 


E(X 4 ) 


( 11 ) 
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Criteria (3), (4), (8) and (10) are applicable to systems with linear damping characteristics. For 
the nonlinearly damped system these four criteria are directly generalized. The conventional 
stochastic linearization technique demands that the mean square difference between the nonlinear 

damping force <p(2f) and its linear counterpart c^X be minimal: 


E 




min 


yielding 

„o> _ £[<P(*)* ] 

" SC* 2 ) 


( 12 ) 


(13) 


In full analogy with the Kazakov’s criterion (5), applicable for the system with nonlinear 
restoring force, we can also require that the mean-square values of the dissipation force and its 
linear equivalent be equal, i.e. 


E [<p ! (i)] = E 


(14) 


which results in 



N 


£ [<p 2 (*)l 

E(X 2 ) 


(15) 


To the best of our knowledge this criterion is formulated here for the first time. The third 
possible criterion of equivalence between the nonlinearly damped system and the one with linear 
damping was formulated by Wang and Zhang (1985) and Elishakoff and Zhang (1992.b). They 

required that the mean-square difference between the energy dissipation function 0(X) of the 

original nonlinear system and that of the equivalent linear system c®X 2 /2 should be minimal 

E [<$(*) - c%X 2 l?§ = min ( 16) 


resulting in the equivalent damping coefficient 



( 17 ) 
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The fourth criterion was suggested by Elishakoff and Zhang (1992.b). It is based on requirement 
of equality of mean squares of <3>(A) and c^X 12, namely 


£[4> 2 W] = 



with attendant value of c M (4) 


c (4) = 2 


N 


E(X) 


(18) 


(19) 


We will evaluate several examples to illustrate the performance of the proposed stochastic 
linearization criteria. 


DUFFING OSCILLATOR UNDER WHITE NOISE 

Consider Duffing oscillator subjected to ideal white noise excitation 

X + $X + aX + e* 3 = W{t) 

The mean-square displacement for the system with e» 0 reads 

E(X%^ = KS/afc el 


( 20 ) 


( 21 ) 


where 5 is the value of the spectral density of excitation W(t). The exact probability density of 
X is uniformly available in the literature and will allow a comparison between the conventional 
and the new stochastic linearization methods: 


Pjx) = C 2 exp 


nS 




( 22 ) 


where C 2 is a normalization constant. In view of Eq.(21), Eq.(22) can be rewritten as 


Pj, x ) = C 2 ex P 




(23) 


We introduce new variables 
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X = pr, 


P = 


( A 2 
4ae 0 


\ 1/4 


\ / 


y = ^(tJ 


,1/2 


(24) 


The normalization condition 


C / exp 


_ 1 

1 , 

1 e A \] 


—X“ + x 'll 

2 

e 0 

l 2 

4 a J 


\dx = 1 


(25) 


yields 


c, = \pzm i 


(26) 


where 


Z.(y ) = 2 /exp( -x 4 - 4yV)ih 
0 


(27) 


Mean square displacement 


E(X*) - !xy/x)dx 


(28) 


reads 

P^OO . „ (29) 

“W •z.wUJ 



where 


Z 2 (y) = 2 /x 2 exp( -x 4 - 4y 2 t 2 )dx 
0 


(30) 


Functions Z t (y) and Z 2 (y) are defined and tabulated for certain values ofy by Stratonovich 
(1961). Note that these functions can also be reduced to cylindrical functions of a fractional 
order (Bolotin; 1984, Piszczek and Niziol, 1984; Constantinou, 1985). Consider a particular case 

e 0 2 = 0.54 , a/e = 1 ( 31 ) 


In these circumstances, 
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y = 0.34021 , ZjO) = 1.26368 , Z 2 (y) = 0.26310 


(32) 


The exact mean square value becomes 

EiX 2 ) = 0.306 ( 33 ) 

Let us now contrast the performance of the conventional and new stochastic linearization 
techniques. Conventional stochastic linearization yields 

kff = a + 3zE{X 2 ) ( 34 ) 

Substitution into expression for the mean-square value 


E(X*) = 3iS/Jfe®P 


(35) 


yields 


or, in view of Eq.(34) 


EiX 2 ) = jtf/p[a + 3 e EiX 7 )] 


E(X *) = e 0 2 /[l + 3(e/a)£(X 2 )] 

resulting in a quadratic 

3(e/a)[£(Z 2 )] 2 + EiX 2 ) - e 0 2 = 0 

For numerical values adopted in Eq.(44), we have, instead of Eq.(38) 

3[, E(X ^j 2 + EiX 2 ) - 0.54 = 0 


with attendant mean-square value 

EiX 2 ) = 0.289 


(36) 


(37) 


(38) 


(39) 


(40) 


which constitutes a difference of 6.47% with the exact solution. 

Consider now the energy based stochastic linearization method. Eq.(6) yields 
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k£> = a + l-SeEiX 2 ) 


(41) 


Substitution k^ 2) instead of k eq <!> in Eq.(35) yields 



E^X 2 ) = JtS/p[a + 2.5e£(X 2 )] 

(42) 

or in view of Eq.(21) 


E(X 2 ) = e 0 2 /[l + 2.5(e/a)£(X 2 )] 

(43) 

This results in a quadratic 


2.5(e/a) [EiX 2 )] 2 + EiX 2 ) - e 0 2 = 0 

(44) 

For the values in Eq.(44), 

we get an equation 



2.5[£(X 2 )] 2 + EiX 2 ) - 0.54 = 0 

(45) 


with attendant mean-square value E(X *) = 0.306 which coincides with the exact value given 
in Eq.(33). For value of e 0 2 in vicinity of 0.54 and ratio e/a in vicinity of unity, the relative 
error furnished by a new stochastic linearization technique may constitute about one percent, 
much smaller than the one produced by the conventional stochastic linearization technique. 

The discovered coincidence of the stochastic linearization result with exact solution 
suggests that for specific set of parameters (namely for e 2 - 0.54) the new version of stochastic 
linearization constitutes a "true" linearization, in terminology of Kozin (1987). 


DUFFING OSCILLATOR UNDER COLORED NOISE 

For the Duffing oscillator subjected to Gaussian white noise the problem is amenable to 
exact solution. Consider now the case where the exact solution is absent. Indeed, the power of 
the stochastic linearization lies in its applicability when all other analytical methods may fail. 
Consider a Duffing oscillator subjected to a colored noise (Falsone and Elishakoff, 1992). The 
motion of the system is governed by the differential Eq.(20), where, under new circumstances, 
Q(t) is a filtered noise; in particular, we assume that Q(t) is the response of the following first 
order filter equation 

Q - - + yW(t) (46) 

W(t) being a white noise with constant spectral density S. In equation (46) the parameter y gives 
the measure of the filtering. In fact, it is easy to verify that the spectral density function S Q ( to) 
of Q(t) is given by 
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(47) 


S Q ( <*>) = 


2 2 

+ co 


and, for large values of y, S Q ( co) tends to S. 

According to the classical linearization technique, Eq.(46) is replaced by the linearized 
one in which the linearized spring k eq (I) is given by Eq.(34). In this way, the stationary mean 
square value response of the linearized system reads 


E\X 2 ] = — * — J tS (48) 

ClH*™ * Pr * V 2 ) 

It is worth noting that, for large values of y, this quantity tends to the mean square response of 
the linear system subjected to the white noise W(t). 

According to the new linearization technique, specified in Eq.(ll), the linearized spring 
Jc (3) is given by Eq. (41). The stationary mean-square response is given in Eq.(48) where k eq 1} 
is replaced by k eq <3) . 

The two linearization techniques have been first applied to the Duffing oscillator by 
varying the filter parameter y and fixing the system parameters at values utilized by (Elishakoff, 
1991): 

jLS/a|3 = 0.54; e/a = 1.00 ^^) 


These values are the same ones for which the new linearization technique yielded the 
exact stationary solution, when the input is a white noise (Elishakoff, 1991). Falsone and 
Elishakoff (1992) evaluated the percentage-wise error between the results obtained by means of 
each of the two linearization schemes, and the results obtained by means of Monte Carlo 
simulation. The superior accuracy of the new linearization approach than that of the classical 
one was evident for each value of the filter parameter y examined. Moreover, it is worth noting 
that, for y = 10, the percentage error of the results obtained by the new technique is practically 
zero, confirming the previous result obtained by Elishakoff (1991). In the paper by Falsone and 
Elishakoff (1992) the percentage error for y = 1, n S/aP = 0.54 and varying the e/a was 
studied. The preferable accuracy of the new approach established even for high level of non- 
linearity; indeed for e/a = 2 the classical linearization yielded an error of 12% whereas the error 
within the new method was under two percent. 

These results confirm that, in the case of a Duffing oscillator, the linearization with 
respect to the potential energy yields much more accurate results than the classical linearization 
approach, even when the input is a colored noise and the level of non-linearity is high. 
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NONLINEARLY DAMPED SYSTEMS 


Let us elucidate now the proposed criteria on an example which was studied earlier in the 
literature via the perturbation method (Khabbaz, 1964). Namely, the system governed by the 
following differential equation 

mX + cJC (1+b^X 2 + + a 4* 4 ) = MO (50) 


is investigated. In Eq. (11) m is the mass, c 0 and k 0 are coefficients modelling the damping and 
spring stiffness of the linear system which is obtained by a formal substitution: a 2 = b 2 = a 4 - 
b 4 = 0. Here coefficients a v b^ a 4 and b 4 are assumed to be positive and specified, W(t) is a 
stationary Gaussian white noise. The exact solution of Eq. (50) is unavailable. We derive 
approximate solutions through various stochastic linearization criteria. The results axe compared 
both with those yielded by the perturbation method and by the Monte Carlo simulation. 

We first consider the analysis through conventional stochastic linearization. Conditions (4) and 
(13) yield respectively, 

= *„[£(**) - a,E (JO * E(X‘)]lE(X*) <51) 


c™ - 4 EiX 1 ) * bMX') - (52) 

The mean square displacement of the linear system (with a 2 = b 2 = a 4 = b 4 = 0 ) equals 

= (53) 

where 5 is the spectral density of W(t). For the equivalent linear system the mean-square 
displacement is 



(54) 


Moreover, the mean-square velocity equals 

<4 = k eq a x /m= ^ /c «i m (55) 

Utilizing the postulated jointly normal probability density of the displacement and velocity for 
evaluation of Eqs. (51) and (52), results in 
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(56) 


= k Q (l+3a 2 o x + 15 a 4 o x ) 

ci, l) = 4l + 3b 2 (kjm)(l + 3a 2 c£ + l5a 4 o 4 x )a 2 x 
+15b 4 (ko/m 2 )(l + 3a 2 a* + 15a 4 c£) 2 a*] 

Substitution of Eqs. (56) and (57) into Eq. (54) yields a polynomial equation 

E A t (cfy‘ = 0 (^8) 


where coefficients A i read 

= -°r,> A i * A 2 = 3 (ajn+tyjlm 

A 3 = 3 (5a 4 m 2 +6mk Q a 2 b 2 +5kQb 4 ) /m 2 

A 4 = 9k 0 (3ma2b 2 +10ma 4 b 2 +15k 0 a 2 b^/m 2 

A s = I35k 0 (2ma 2 a 4 b 2 +3k 0 afb 4 +5k 0 a 4 bj/m 2 

A 6 = \35k 0 {5ma^b 2 +3k Q a2b 4 +30k 0 a 2 a 4 by) /m 2 

A 7 = 2025k 2 a 4 b 4 (3aZ+5aJ/m 2 

A s = 30375 k^a 2 a^b 4 jm 2 , A, = 50625 ^W" 1 2 


Here ay is the mean-square response of the corresponding linear system. According to the 

Descartes’ rule of signs, the polynomial Eq.(58) has a single positive root for o x , since the 
sequence of coefficients in Eq.(18) has only one change in sign. The evaluation of the mean- 
square displacement should be performed numerically. 

Let us now utilize the criteria (8) and (16). The potential energy of the system reads 

U(X) = k 0 {(II2)X 2 + (1/4 )af* + (1/6) aJC 6 ) (60) 

For the energy dissipation function we obtain 
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( 61 ) 


<D(A!) = / c 0 u( 1 + b 2 u 2 + b A u*)du 
o 

- c„{(l/2)Jf' * (1/4)6/' * (1/6)6/“} 


The expression for k xq (3) reads 


- 2*o K, 


( 2 ) 


where 


- (1/2) * (S/4)«,o 5 - (35/6) a 4 o' 
The expression for the equivalent linear damping reads 


(62) 


(63) 


c«> = 2 c„ [(1/2) * (574)6,0* ♦ (35/6)6 4 o'J - 2cf, 


( 2 ) 


where 


-(2) _ 


1 5, 2^/i 5 2 35 4 \ 2 

— + -i> 2 — + -a 2 o x + —a o x o x 

2 4 m 2 4 6 


35 , 4£ 0 / i 5 2 

\ — + — a.o x 


35 


m 


2 4 


4) 2 4 

VM °* 


(64) 


(65) 


The mean-square displacement equals 


CT r = 


jlS 


a„ 




< 2) < } 


( 66 ) 


Substitution of Eqs.(63) and (64) into Eq.(66) leaves us with a polynomial equation 

i 5<a 2 )‘ = 0 


(67) 


where the coefficients B, are defined as 
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B l = 1, B 2 = 5{ajn + IcJbJ/lm 


B o = -<V’ 

B 3 = 5( 14a 4 m 2 + 1 5mk 0 ap 2 + 14k^b^6m 2 

Z? 4 = 25 £ 0 (15a 2 2 £> 2 m + 56a p 2 m + 84k 0 apJ /24m 2 

= 175 fc 0 (10 a 2 a 4 b 2 m + 15£ 0 a 2 2 £> 4 + 28k 0 apJ/12m 2 ( 68 ) 

5 6 = 875 fc 0 (2$a A b 2 m + 15* 0 a 2 3 /> 4 + I68k 0 a 2 apyi2m 2 
B n = 6125 ^ 2 a 4 f» 4 (15a 2 2 + 28a /36m 2 
fl 8 = 214375 kja 2 a?b 4 /18m 2 , B g = 1500625 * 0 2 a 4 3 Z? 4 /81m 2 

Again, Eq.(67), like Eq.(58) has a single positive root. Combined utilization of criteria (10) and 
(18) yields 

C - 2**? < 69 > 

where 


Also 


where 


< 3) = 


-1 + — a,ov + 35 
4 4 2 * 


/ 1 


105 


16 2 6 
. 11/2 


W* 


105a 4 o x 


C (3) = 2c C (3) 

^eq z - L o'-cq 





r I 2 

1 5 

+ 



A: (3) 

4 4 

m x (16 2 

6 J 

m 


105 , , 

+ -T- b 2 b * 


fc® 


T3 


m 


a x + 1056 4 


t (3) 

*7 


m 


a. 


1/2 


(70) 


(71) 


(72) 


Substitution of expressions for equivalent spring stiffness and equivalent damping coefficient into 
Eq.(54) yields an equation for o x 2 : 
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(73) 


2 

Ox 


xS 

-Wi-P) 

L eq K eq 


2 

qr 


2 

< 3 ) ci 3) 


= 0 


The explicit form of the resulting polynomial equation in terms of a x 2 is cumbersome and is not 
reproduced here. One can show, however, that the resulting equation also possesses a single root, 
as in previous cases. 

The numerical results for Eq.(58)~(67) have been obtained by Elishakoff and Zhang (1992.b). 
Comparison with the results of Monte Carlo method (see Fig. 1 in the above paper) demonstrate 
that the conventional stochastic linearization technique results in the largest error. For small 
values of the parameter e criteria 3 and 4 yielded values which were extremely close to the 
simulated meau-square responses. For intermediate values of e, namely e - 1, the fourth criterion 
performed the best, whereas for larger values, namely for e - 2, the third criterion yielded resu ts 
in close vicinity to those of the Monte Carlo method. Elishakoff and Zhang (1992.b) considered 
also a softening system: 


mX + cX + (k 0 /e)sgn(X)(l - e^) = W(t) 


For values of e up to unity the criterion of equal energy variances yielded results in best 
agreement with the simulation. In the range Ues2 the energy wise minimum mean square 
criterion turned out to be superior to other criteria. 


HYBRID WEIGHTED STOCHASTIC LINEARIZATION-MONTE CARLO METHOD 


The motivating considerations for developing a hybrid stochastic linearization-Monte 
Carlo Method are as follows. 


The fraction of problems amenable to exact solutions is very small. For most problems 
the exact solutions are unavailable. In these circumstances either purely numerical approximate 
techniques are utilized, or the Monte Carlo method is applied. Amount of computations within 
the Monte Carlo solution of the problem may be enormous for a large system. One still may want 
some analytical method combined with small scale simulation. One can significantly reduce the 
amount of calculations by the proposed method. To do this we choose the following form o 
the weighting functions (Elishakoff and Colombi, 1993) 


H<*) 


/ 


l 1 


xi/2 


+ laU 1 
*'■ 1 ‘ / 


(75) 


where coefficients and p, should be determined from the numerical experiments by the Monte 
Carlo method; n signifies the number of series of Monte Carlo simulations. Note that weighting 
functions have also been considered, although in totally different contexts by Wang and Zhang 
(1985), Yu and Cao (1988), Izumi et al (1989), Fang and Fang (1991), Elishakoff and Zhang 
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(1991.a) and Zhang (1992). 

In other respects however, the proposed method is that of stochastic linearization. The four 
alternative criteria are replaced by their equivalents utilizing the weighting functions. For 
example, the third, criterion is replaced by 

e[mX)\jJ{X) - k^X 2 /. 2 ]} = min (?6) 

yielding 

_ 2£MaW(A)* ; 1 (77) 

"* E{w\X)X 4 } 


Each simulation series should be conducted for specified sets of parameters of the system. 
Thus for n specified sets of parameters the results of Monte Carlo simulations will numerically 
coincide with the results of the stochastic linearization technique. It is expected that for other 
sets of parameters the accuracy of the results yielded by the stochastic linearization technique 
with weighting function will be quite satisfactory. 

Consider, for illustration purposes, a system governed by the following differential 
equation 

mX + cX + kjX + kJC s = W(t) ( 78 > 

where m, c, k 2 and k 2 are positive constants, W(t) is an ideal white noise with zero mean. The 
system (82) is amenable of exact solution. It has been chosen in order to elucidate the errors 
associated with approximate techniques. It was shown by Zhang, Elishakoff and Zhang (1991) 
that energy based linearization technique, even without recourse to the weighting functions 
reduces the errors in determining the mean-square displacement of the system, by about 50%. 

We will illustrate the application of the energy based stochastic linearization technique 
with weighing function. The potential energy of the nonlinear system in Eq.(83) reads 

U(X) = (1/2 )k r X 2 + (V6)k^ 6 (79) 


We perform a single series of Monte Carlo simulations, 
at unity, with 

MX) = (1 + a 1 L/) 1/2 


Hence, in Eq.(79), n is fixed 

(80) 


Substitution of Eq.(112) in Eq.(77), in view of Eq.(lll) results in 

El 6X S + 3a k.X 10 + a k^ 14 ] 

Jc ~ Ic + 2k . 

1 2 £[363 / 4 + 18c tkj 6 h- 6a kJC 10 ] 


(81) 
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The spectral analysis, yields the following mean-square displacement 


2 JiS °Jr. 

Ov = ■' “ — 

ck 1 +A 

“t 


where ot is defined as 


OX' - ns /ck l 


The "corrective" term A reads 

^ k 2 E[6X* + 3o kjt 10 + a kjt 14 ] 

T l ~E^6X r Tl2oJc^ r ^Ja^ 


( 82 ) 


(83) 


(84) 


To calculate E[X 2 ’] we use the usual approximation that X has a normal probability 
density with zero mean Eq.(88) becomes: 

15015o ik}o x + 630o* 2 c£ + {l0k 2 - 'il5akjc l a x ^o x ^ 

+ ISaJcyOx 46^ - 15ak 2 a x o x )c£ - 6 k^ = 0 


In the latter equation the parameter a is not known. In order to determine it, we first 

solve an auxiliary problem of evaluation the mean square displacement o x for specified set of 
parameters, by the Monte Carlo method, say for 


m 


= m w 


- 


= c<», k, = fc"' 


( 86 ) 


The result of simulation is denoted by a x . This allows one to determine the value ofa = a 
corresponding the simulation results 

- li^Skfix - 3 k y o x + 3^,0,^ 

15 Oy( 1 00 1 ^ a x + - 21^6^ + k?o x - k 2 a x a x ) 

We fix a at the value a determined from this equation, and use it in Eq.(85) for values 
of parameters, other than those listed in Eq.(86). Once dt is substituted, Eq.(85) becomes a 
polynomial equation with respect to a x : 
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35035k? a]? - 35035k? a™ + 4473/t 1 /t 2 2 a“ + 3003 kJc?o x 

- 1470 k?o x ~ 3738 k ] k?o x o x - 2268 k ] k?o x o x - 35k?k 2 o x o 7 x 

+ 161 k?k 2 o 7 x + 35k?k 2 o x o x - 35k?k 2 o x + 126 k v k 2 o 6 x - 189 k?k 2 o x o 5 x 

’ ° ( 88 ) 

- 63 k?k 2 o x o x - 126k l k 2 a x a x + 63 k?k 2 a?o x - 3k?o x o x 
+ 3k? o? + 63k{k 2 o 2 x p 2 x - 3 k?o x o x + 3fc t 3 a* + 3k?o 2 x p x 
- 3k?o x o x + 3 k?o 2 X ' - 3 k?o x - 0 


To gain some insight, let us consider some numerical results. Let 
m (1 > = 1 , c (1) =0.1 , k? l) = 10 , & 2 (1) = 15. Simulation results are in close vicinity with exact 

solution. Calculations yield o x - 0.687. The conventional linearization yields for this set of 
parameters o x = 0.5043, or 23.44% off the exact solution, and about 20% off the simulation for 
a sample of 10 6 simulating systems. The energy based stochastic linearization without weighting 
function results in an estimate o x = 0.5456, or 17% off the exact value. 


The value of &, matching the results of the Monte Carlo method and stochastic 
linearization equals &= - 0.000607. As noted above, at k^—10 the stochastic linearization 
technique yields results coincident with those furnished by the Monte Carlo Method. At value 
k? 2> = 11, the conventional stochastic linearization is off the solution predicted by the Monte 
Carlo method by 22.77%. The energy based criterion without weighting function results in 17% 
of error. The proposed combination of the energy-wise linearization with the Monte Carlo 
method yields an error under 5%, which corresponds to the reduction of error associated with 
the conventional stochastic linearization, by the factor of 4.6. In the closer vicinity to the 
parameters for which Monte Carlo analysis was performed error may be significantly smaller. 
For example, for k/ 3> = 10.5, the proposed method yields an error of only 3.85%. 


CONCLUSION 

The results of the application have confirmed a superior performance of the new linearization 
techniques, for various cases, including both "hard" and "soft" oscillators. It appears that 
additional studies will prove useful to study the applicability of the proposed new stochastic 
linearization criteria in civil, mechanical, aeronautical, aerospace and naval engineering contexts. 
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Abstract— A new stochastic linearization technique is employed to investigate the large amplitude random vibrations of 
a simply-supported or a clamped beam on elastic foundation under a stochastic loading which is space-wise either (a) 
white-noise or (b) uniformly distributed load and time-wise white noise. The new version of the stochastic linearization 
method is based on the requirement that mean square deviation of the strain energy of the nonlinearly deformed beam, 
and the strain energy of the equivalent beam in a linear state, should be minimal. As a result, the modal mean square 
displacements are expressed as solutions of a set of nonlinear algebraic equations. Results obtained by conventional 
equivalent linearization method and by the new technique are compared with the numerical results obtained from 
integration of the exact probability density function (when exact solution is available) or with the result of the Monte 
Carlo simulations (when the exact solution is unavailable). It is shown that the new stochastic linearization technique 
yields much more accurate estimate of the mean square displacement than the classical linearization method, which has 
attracted in the past interest of about 400 investigators in variety of nonlinear random vibration problems. 
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INTRODUCTION 


Since it was first proposed by Booton [1] and Kazakov [2] four decades ago, the 
stochastic linearization method has been applied in numerous engineering fields. Application of 
this method to mechanical systems was performed by Caughey [3-5]. The excellent reviews of 
these studies have been given by Roberts and Spanos [6] and Socha and Soong [7]. Several 
authors recently have studied the various new versions of the improved stochastic linearization 
technique [8-10], which can be classified into two main categories: energy-based versions [8] and 
weighting function based versions [9-10]. Review of these developments was recently presented 
by Elishakoff and Falsone [11]. In the recent study [12] by the present authors, the beam under 
the time-wise and space-wise white noise loading was studied by the conventional and the new 
version of stochastic linearization technique, which is based on the criterion of minimum mean 
square deviation of potential energies. Comparison was performed with the exact solution 
available in that case. In the present study, the same energy-based method is employed to 
investigate the beam under general time dependent stationary random excitation, when exact 
solution is unavailable. The beam is resting on an elastic foundation and simply supported at 
both ends, and loaded by a stationary random excitation. As a result of moderately large 
vibrations, the different modes are coupled with each other, whose effect has been shown quite 
significant in Ref. [13-14]. Numerical simulations are carried out to compare their results with 
those yielded by the new and conventional linearization techniques. It is found that the new 
version of the stochastic linearization technique yields considerably more accurate results for the 
mean square response of the beam than the conventional stochastic linearization technique, 
especially in the case of large nonlinearity. 
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PROBLEM FORMULATION 


Consider a beam resting on an elastic foundation, restrained at its ends, and loaded with 
spacewise distributed time dependent loading p(x,t). The equation governing moderately large 
vibrations reads [8] 


( 1 ) 


c .,d 4 w .rtf 2 ™ . dhv a dw „ , A 

El - N + p A + 6 + K.w = p( x, t ) 

dx 4 dx 2 dt 2 dt f 


where the axial force is given by 


N = 



( 2 ) 


and E is the elastic modulus, /= moment of inertia of the cross section, p=mass density, A=cross 
sectional area, Z,=length of beam, ^=stiffness of the elastic foundation, p(x,t )= time dependent 
stationary and ergodic random distributed loading acting on the beam. Moreover, for the load 
we will make an assumption of the multiplicative representation, i.e. 


p(x, t ) = r(x)q(t) 


( 3 ) 


where r(x) is the deterministic function of the axial coordinate only and q(t) is the random 
function of the time alone. Also the autocorrelation function of the loading, which is considered 
to be weakly stationary in time, is known. According to the Wiener-Khintchine relationship 


- >i) - 


where S(oi) is the spectral function of q(t). 
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Let the beam deflection function be expanded in terms of appropriate orthonormal set of 


mode shapes in vacuum § n (x) of the associated linear structure as 


w 


= E W M ♦.(*) 




( 5 ) 


where ty n (x) satisfy the following relationships 


d\ j> 

El __1 = pAco„$ 
dx 

1 

= 6^ i = x/l 

Then the expression for the axial force becomes 


( 6 ) 


N = 


2L m " dx dx 


( 7 ) 


EA 

2L 1 


YYd ww 

^ mn m n 


m*l n»l 


where 


D = 




l 


d% d% 




( 8 ) 


In view of Eq.(6), the governing equation (1) takes the following form 
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E 



JL W + co’W + — f -W)<k n 

pA " H pA " 


-YY d ww.w 
4 Z-* i] 1 J n 


2p L i*l ;*1 


<^ 2 


= r(x)^(r) ( 9) 
pA 


Multiplication of Eq.(9) by <j> m , integration over the length of the beam and the use of 
orthogonality conditions given in Eq.(6) yield a set of coupled nonlinear differential equations 
for the modal amplitudes WJt) 


w * JLw. * «>„ * 5lw. * * -tiio^wv/w - ^K<) ( 10 ) 

P A 2p L ,*1 1 rt-1 


£ 

pA 


where 


i 

c.=4*.®k i)^ 



cm 


By either conventional or the new version of the stochastic linearization technique, Eq.(10) is 
replaced by a linear one 


W 


JLw 

A m 

pA 


+ k o i 2 m W m 

m,eq m rn 



(12) 


where k is the nondimensional equivalent stiffness. The problem consists in evaluating k meq 

ntyCq 

through two versions of stochastic linearization technique and comparing the results for response 
with those obtained from the Monte Carlo simulations. 
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ANALYSIS 


The potential energy of the system, represented by Eq.(10) reads, in conjunction with Eq. 


( 9 ) 


m oo 00 00 


U = jr I(a.i - %vl * E- 4r \(L'LY, D >i R ~ w : w , w .') dw - (13) 

m-l 2 P A m*l 2pL ‘" l V* 1 


The main idea of the modified method consists in the requirement that the mean square 
difference between the potential energy, associated with the original nonlinear equation, Eq.(10), 
and its equivalent linear counterpart Eq.(12), to be minimal. That is 


U -Y -k a i„W, 

n,eq n ' 

n-L ^ 


= min 


which is achieved by demanding 





r 

2 

£ 

“■Ek,* 
««1 ^ 



= 0 


O = 1,2,-) 


The conditions (15) are equivalent to 


(14) 


( 15 ) 


U - Y —k Oil W 

/ -j ~ n,eq n r 


w: 


= 0 


(m = 1, 2, -) 

Thus, if only first N normal modes are included, Eq.(16) can be written as 


(16) 
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( 17 ) 


[A](y = 2iB) 


where 



WX 2 ] 


... E[W 2 W 2 ] 

[A] = 


£[w> 2 2 ] 

... 



E[wZw 2 2 ] 

... E[W^Wi\ 

{k 7 = ik. 
« ? 1 . 

Ic "• 

• W 


(5F = {E[W?U] E[W?U] - E[W£U]) 

The solution can be written as 

k Ntq <3) N Y = 2 [A]' 1 {B) 

For the Gaussian random processes with zero mean, we have [14] 


£[ Z 1 Z 2- Z J = E [II £ [¥J 

a// independent pairs > }+k 


(18) 


(19) 


( 20 ) 


where the number of independent pairs is equal to (2m)\/2 m m\. In view of the representation of 
the potential energy by Eq.(13), the following terms will appear in {B} 
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E[W 2 m Wt] =E[W 2 m \E[W 2 ] * 2(E[WW n ]) 2 
E[WfW 2 m Wl] = E[W 2 ]E[W 2 W 2 ] h- 2E[W^ m )E[W,WW^ 

* 2E[W l W n ]E[W,W 2 W n ] ( 21 ) 

= E[W?]E[W 2 W 2 } + 2(£[W,WJ) 2 £[W B 2 ] 

+ 2(£[W,WJ) 2 £[W m 2 ] + 8£[W,WJ E[WyV n \ E[WW n \ 

The moments E[W m W n ] are related with the spectra of the system in the following way 

OO 00 

E[ww n ] = J f hjt -g *,(* -g ^wo^gi 

\P*^/ -00-00 


( 22 ) 


where the function hjt) is the response of the SDOF system represented by each of equations 
(12) to a unit impulse. The Fourier transform of hjt) 


oo 

AJ 0 - 2- J 

•OO 

00 

- J e hj, t)dt 


(23) 


is given by 


«■» * 


I 2 ? 

k - or 

m, « 




co 


pA 


For stationary random excitations 


(24) 


E[q,(t)q 2 (t)} =R(t 2 ~t Y ) = 


( 25 ) 


7-8 


With the use of Eq. (23), the double integral in Eq. (22) becomes 


= ff h S t ~K) bnit-Q d h dt i 

- 00-00 


- jH m (-u)HJo>)S(o>)dm 


( 26 ) 


For the white noise excitation, can be evaluated through frequency response function of 
system (12) by means of residue theorem to yield 




5 = 


pA 


(k at 2 -k co l) 2 + 2(JL.) 2 (k o o 2 + k to!!) 

V «, eq ** n, eq n ' V pA. ' ' m - m n. ea n ' 


(27) 


n, eq 


m, eq 


n , eq 


where S 0 is a constant value of the spectral density S( oj). As a result, is expressed through 
the sought parameters k^ eq . After substituting Eq.(27) to Eq.(22), and the result into Eqs.(21), 
(19), (18), we obtain a set of nonlinear algebraic equations, i.e, Eq.(17) in terms of k , which 
must be solved numerically to obtain k . 

From the conventional stochastic linearization technique, we have 


k = 1 + 

m,eq 


K, 




E (^-) ! E £ E d,«„Q,Q,Q. <S ‘ ,S ~ 1 2S -^ (28) 


P Ao£ 2p L'vte. P 


The overbar denotes the association with the conventional stochastic linearization. Note that in 
obtaining the latter expression for k m , we utilized the same procedure as Seide [15]; indeed 

when Kf= 0, Eq.(28) coincides with Eq.(22) of Ref. [15]. In the same fashion as for the new 
stochastic linearization, we have arrived at a set of nonlinear algebraic equations which must be 
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solved numerically. 

A particular case of the beam simply supported at both ends under a spacewise uniform 
random excitation is numerically evaluated: 


<L(D = t/2* sin /mi| 


Q m = 


D = 


R = D 


0 

2\j2 

rrm 

0 


El 


if m is even 
if m is odd 

if m = n 


m 2 7t 2 if m - n 


2 X-r 1 4 4 

o)_ = mrr 


pAL 4 


Then accordance with Eq.(13) the potential energy reads 


(29) 


where 


O3 0 

U = _ 

2 


* °E «'„ 2 * -2 






4i? 


£ "* 2 ^ 2 

, m-l j 


a:. 


a = 


PAcOq 


tOn = 


n 4 EI 
p AL 4 


2 2 4 
= to 0 m , 


R = 


N 


~A 


(30) 


(31) 
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Therefore Eq. (19) becomes 


K* »«**., Ar ‘^.^} r - 4 W‘ { £ I "VH } r (32) 

O ) 0 

This set represents a set of nonlinear algebraic equations since E[W^V\ themselves contain k. . 

For the simply supported beam with elastic foundation, the motion of equation (10) 
becomes 


W + 


JL W + 0) 2 m 

p A m 


1 + 


a 


m 


4 2 R 2 m 2 


E '* 2 ^ 2 


IT = 

pA 


(33) 


As a result, the equivalent coefficients in Eq.(28) by the conventional stochastic linearization 
technique become 


k = 1 + 

m, eq 


a 




Ql (5 5 + 25,1) 

\ mm n/i *** ' 


(34) 


m 2R 2 m 2 n .i (pA) 


For R -* oo the contribution of nonlinear terms disappears. Therefore, R in actuality is a 

nonlinearity parameter. We will be interested in evaluating the effect of this parameter on the 
mean square response. 

NUMERICAL RESULTS AND DISCUSSION 
Calculations are performed to compute the mean square responses by the new and 
conventional stochastic linearization techniques. Since no exact solution is available for the beam 
under spacewise excitation, the numerical simulation is also carried out for all the cases of 
parameters listed at the figures. 

In carrying out the numerical computations, the following material parameters have been 
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assumed through out the present work, oo 0 2 =1.0, and P 0 =P/(pA)=0.1. Three normal modes are 
considered in the numerical evaluation, but only the first mode response is plotted in the figures 
since the first mode response dominates the mean square displacements of the beam. In Fig. 1, 
the mean square responses by conventional and new linearization and Monte Carlo simulation 
versus the nonlinearity coefficient R are shown for the case of S o =1.0, a=0.0, (1=0.1. Fig. 2 also 
shows the three results obtained both through conventional and new stochastic linearization 
tec hni ques, and Monte Carlo simulation for another set of parameters with 50=5.0, a=1.0, (1=0.1. 
It can be seen that the new version of stochastic linearization method exhibits a superior 
performance over the conventional one at all ranges of variation of R. Fig. 3 displays the 
influence of the stiffness of the elastic foundation on the mean square response of the system for 
the case 5„=0.1, 0=0.1, R= 1.0. These results clearly indicate that the present method yields 
results in closer vicinity with the simulation results than the conventional linearization technique. 
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Fig. 1 The mean square response of the first mode E[W t 2 ] calculated in the three mode 
approximation (N=3) vs. the nonlinearity coefficient R (S 0 — 1, ct— 0, P—0.1). 


Fig. 2 The mean square response of the first mode EfW^] calculated in the three mode 
approximation (N=3) vs. the nonlinearity coefficient R (S 0 =l, a=l, p=0.1). 


Fig. 3 The mean square response of the first mode E[W t 2 ] calculated in the three mode 
approximation (N=3) vs. the nonlinearity coefficient R (S 0 =5, ct— 1, P—0.1). 


Fig. 4 The influence of the stiffness of the elastic foundation a on the mean square response of 
the system (S 0 =0.1, p=0.1, R=l). 
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Abstract: The new stochastic linearization technique is employed to investigate the nonlinear mean square response of 
a beam with arbitrary boundary conditions under time dependent stationary random excitation. To demonstrate the 
accuracy of this new method, mean square response is obtained through both the new and conventional versions of the 
stochastic linearization technique. An example of a beam with both ends simply supported, is presented, for the case of 
white noise excitation, to illustrate the present method. Numerical simulations are performed to check the accuracy of 
the modified technique. It is shown that the modified version can yield more accurate results for the mean square 
response of the beam than the conventional stochastic linearization technique, especially when there exists a large 

nonlinearity. 

Keywords: conventional and new stochastic linearization, nonlinear beam, mean square response, 
Monte-Carlo simulation 
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1. INTRODUCTION 


The random vibrations of beams in linear and nonlinear settings have been investigated 
by several authors. Linear random vibrations have been investigated by Eringen (1957), 
Bogdanoff and Goldberg (1960), Crandall and Yildiz (1962), Elishakoff and Livshits (1984) and 
Elishakoff (1987), by employing the normal mode method. Eringen (1957), Elishakoff (1987) 
and Elishakoff and Livshits (1984) were able to sum up the infinite series of modal contributions 
and derive closed-form solutions for simply supported beams subjected to loading which is both 
time-wise and space- wise white noise. For this particular excitation, Herbert (1964, 1965) 
succeeded to obtain an exact, although not a closed-form solution, for probability density function 
of modal displacements. For the general case of excitation, when dealing with the nonlinear 
stochastic boundary-value problems, most investigators have employed approximate techniques: 
either the classical perturbation method or the stochastic linearization technique. The latter 
method has attracted numerous investigators. Indeed the only monograph in this subject, that by 
Roberts and Spanos (1990) lists approximately 250 papers utilizing the stochastic linearization 
technique. The review by Sinitsin (1974) lists about 120 studies predominantly performed in 
Russia. The review paper by Socha and Soong (1991) lists numerous publications written in both 
West and East. Thus presently there are approximately 400 studies based on the classical 
stochastic linearization technique. Unfortunately, these methods have severe limitations: the 
perturbation method usually only leads to results of acceptable accuracy for the case of small 
nonlinearity; stochastic linearization technique can yield an error as large as 20% in the estimate 
of the mean-square response for some cases of nonlinearity and excitation. It should be bom in 
mind that even an exact solution by the Fokker-Planck equation method may involve a large 
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amount of multiple integrations to evaluate the mean square responses if many modes need to 
be included. The latter difficulty is of purely numerical nature whereas the approximate methods 
have their inherent difficulties. Furthermore the fact that the (numerically cumbersome) exact 
solution is available for extremely specific cases of excitations rules out its general application. 
The above disadvantages of the approximate methods and absence of the exact solution for the 
general loading case encouraged the present authors to seek for alternative approximate technique. 
Improved stochastic linearization method seems to be an attractive method in this respect due to 
the fact that it not only retains the advantages of the conventional stochastic linearization method 
such as simplicity and straightforward manner of derivations, but also may greatly improve the 
accuracy. Several authors recently studied the various versions of the improved stochastic 
linearization technique (Elishakoff and Zhang, 1991; Elishakoff 1991; Zhang, Elishakoff and 
Zhang 1990; Fang and Fang, 1991). In this study a new stochastic linearization technique, as 
discussed in several references (Elishakoff and Zhang, 1991; Elishakoff, 1991; Zhang, Elishakoff 
and Zhang, 1990), is extended to treat random vibrations of the nonlinearly deformed beam. The 
main idea of the new method consists in the requirement that the mean square value of the 
difference of potential energies of deformation, associated with the original nonlinear equation 
and its equivalent linear counterpart, should be minimal. It is instructive to first elucidate the 
basic idea on the example of the single-degree-of-freedom system, governed by the following 
differential equation 

mX +cx +g(x) = F(t) (*) 

where F(t) is a random excitation resulting in stochastic response x(t); g(x) is a nonlinear 
deterministic function of displacement x. Within the stochastic linearization scheme, this 
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differential equation is replaced by the "equivalent" linear equation 

mx + cx + k^x = F( t) (2) 

where the coefficient is determined through some suitable criterion of equivalence. In the 
linearization scheme utilized by Elishakoff and Zhang (1991), the equivalence criterion is chosen 
as follows 


E 


[U{x)-^k^f 


L = min 


(3) 


where U(x) is the potential energy of deformation of the original nonlinear structure, i.e. 


U(x) = J r Jg(a)da 


(4) 


This is accomplished by requiring 



[U(x)-±k x'Y 


*9 


= 0 


Eq.(5) results in the following expression for the equivalent spring stiffness 


(5) 


k _ 2 E[xHJ(pcy\ (6) 

** E(x *) 

In recent studies (Elishakoff and Zhang, 1991; Elishakoff, 1991; Zhang, Elishakoff and Zhang, 
1990), the authors have demonstrated the accuracy of this linearization technique by comparing 
the computed mean square displacements from different stochastic linearization methods with 
some known exact solutions. In present study, the authors extend the above technique to the 
continuous structures. The main idea is same as the one described for the single-degree-of- 
freedom system, except that the original continuous nonlinear stmcture is replaced by a multi- 
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degree-of-freedom linear system, and a set of equivalent spring stiffnesses are expressed by 
equations analogous to Eq.(6) with x now corresponding to different modal displacements. The 
procedure will be elucidated in detail for random vibrations of the nonlinearly deformed beam. 

In this paper we consider beams simply supported or damped at their ends. Two loading 
conditions are considered: (a) the space-wise and time-wise white noise, in which case the exact 
solution is also obtained, (b) the space-wise uniformly distributed load and time-wise white noise, 
in which case no exact solution is available and the Monte Carlo simulations should be 
performed. In all cases and wide variety of levels of excitation the proposed method turns out 
to be superior to the classical stochastic linearization method. 


2. FORMULATION OF THE PROBLEM 

Consider a beam on elastic foundation with pin-ended supports that are restrained from 
axial motion (Fig. 1). 

The beam is under a loading q(x,t) which is space-wise and time-wise white noise with 
the following auto-correlation function 

R J? C V X 2> 0 = 2jLS 0 6 ( X 2 -*l ) & ( f 2 ~ 0 C 7 ) 

The deflection is represented by the Fourier series in terms of mode shapes of the undamped 
beam 

w(x,t) = £ W n (t ) ( 8 ) 

««i *•> 

w n (t) is the modal contribution corresponding to nth mode. It is assumed that only the first N 


8-5 



modes of the beam are significantly contributing to formulating the response. However, it should 
be bom in mind that the assumption that the power spectral density of the load is that of white 
noise implies that all the modes are excited and contribute to the response of the beam; N is 
determined by the required accuracy in evaluation of the specific response characteristic, such 
as mean-square displacement or mean-square stress. Crandall and Yildiz (1962) have shown that 
for the linearly deformed beam under white noise excitation, if the infinite series representing 
quantities such as mean square displacement, mean square stress, etc., converge, then the results 
can be made as accurate as desired by taking sufficiently large N. It is reasonable to expect 
similar results for the nonlinear problem. We consider a Bemoulli-Euler beam with transverse 
damping. Due to the fact that the equivalence criterion will utilize the concept of energy, we 
first formulate appropriate energies of the beam. Kinetic energy of the beam is given by, with 
Eq.(8) taken into account 


T = dx = 

2 Jo { dt) 4 


the potential energy of bending deformation reads 


( 9 ) 


EI rl 

, \ 2 
dhv 

T Jo 

l 3 * 2 ) 


. n*EI 4* 
dx = 22 n 


AL 3 


4 2 


( 10 ) 


Potential energy of stretching is given by 


V - 


AE 

~2L 


ir( ^dx 

2Jo (dx 


7l 4 AE 
"32 U 




if 


S- 

/l«l 


2 2 


(ii) 


whereas the potential energy of the deformation due to elastic foundation is 
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v . = l" j k f w2dx = ■jV'E w - (12) 

where k f is the translational stiffness of Winkler foundation. Potential function of the load is 


v i = ~J r Q L q(x,t)iv(x,t)dx 


We expand the load in the series analogous to Eq.(8) 


N 


MIX 


<fat) = £ ^( f )sin_ 

where again N terms have been retained. Then V) becomes 


(13) 


(14) 


N 


v i = -yE 

The Lagrangian ££=T-V, where V- V b +V s + V e + V b may now be written as 


(15) 


g = (wf 

4 n) 4 L 3 £( 


4 „.i 
x 4 EA 


32L 


N 


12 


s» 

w*l 


2 2 
w« 




if 




If 

t e <?, 


w 

h m 


The equations of motion are 


(16) 


d_ 

~dt 


' d<£ N 


dvv 

V "/ 


dg 

dw 


= 0 (n = 1,2, -,1V) 


(17) 


where w, are considered as generalized coordinates. Substitution of Eq.(16) into (17) yields the 
following N equations in terms of the modal amplitudes w m (t) 
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Eq.(18) is a nonlinear stochastic differential equation. We seek the mean square response of the 
modal amplitudes. In this study, a new stochastic linearization technique described in the 
preceding section for the single-degree-of-freedom system, is generalized to investigate the 
continuous structure at hand. The nonlinear system (18) is replaced by the following equivalent 
linear one 

W n +-£-vv„ = fjx) (n = 1, 2, -,N) (20) 

pA 

In Eq.(20), it is assumed that the replacing linear system is a decoupled one. The question of 
the decoupling will be addressed in detail in Appendix. 

In our problem, the total potential energy of the system (18) is 
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u ( w v *" w n) = 


n*EI 


N 


N 


2pAL 




<4 A 

- t E" 
L *- 1 


T/i\ 2 +_2_V 
4 • 2 P a^t 

2 AT 

-E 


w 2 + 


ji 4 £ 


12 


4 vv 2 


2 

a<O n JL 2 

W„ + 


2 

O) 0 




16/? 


16 pi 4 

E" 


E« 


V 




( 21 ) 


12 


2 w 2 


#i-l 


We generalize the requirement of Eq.(3) valid for the single-degree-of-freedom system, for 
continuous beam by requiring 


12 


AT 


U(w v w 2 , •••, w N ) - £ w 2 




= min 


which is achieved by using conditions 


( 22 ) 


d \ 
dk i m) 


W 2» ,"%MV) - 53 

»-l ^ 


= 0 


(23) 


(m = 1, 2, JV) 


The conditions (23) are equivalent to 







N 1 


E' 


U(W V W 2 , ,-,w N ) - 53 ±k£w? 

2 


= 0 


(m = 1,2, -,N) 


(24) 


After some algebra, we arrive at the following expression for nonlinear spring constants 


(*„} =2[A]-‘W 


(25) 


where 
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U„F = {*<■’ k% - C} 

(SF - {£[w, ! U] E[w^U] £[hvC7]} 


We denote 


E[h' 1 2 *' 1 2 ] £[h' 1 2 w 2 2 ] - 

£KV] £[w> 2 2 ] - 


[A] = 


£[w 2 2 w 2 ] 


( 26 ) 


pK^w 2 ] £[^w 2 2 ] - £[*vh$ 


y. = £[w 2 ] ( 27 > 

Gaussian assumption for distribution of vv ( - and the subsequent conclusion that the equivalent 
system is decoupled (see Appendix), leads to the independence between different modal 
amplitudes w ; and vv ; (i*j). Therefore 


E[w?) = 3y 2 

£[w ; 2 w ; 2 ] = y.y ; . i*j 

For simplicity, let us investigate the particular case N=3. We have 


[A]’ 1 


1 

10 >W3 


* 

-y 3 

—y 2 


4 w 3 

yi 

-y i 

-y a 

-y i 

>3 


(28) 


(29) 
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Consequently we obtain 


£K ! U] - «*»-a * * -^T 

2 -.1 2 m mi 10 K 


(30) 


Under the assumption that the system is driven by zero mean Gaussian white noise q„ with 
spectral density S„, i.e. f„ with spectral density SJ(pA) 2 , we obtain from Eq.(20) 

t<-> - 4*1/ y. < 31 > 


<*„l * oJu<;( l/y, 1/y, - 1 fy n } T 


where 


2. L*S 0 
°° 


Substituting kj n> in Eq.(32) into Eq.(25) and noting the fact that w H is Gaussian yields 


to! UOW: 


4y 2 y 3 -yji ~yj 2 £ K 2 ^ 

-y 2 y 3 *yj 3 ~y^2 ■ e[w 2 u] 

-y 2 y 3 -yj* e[w?u] 


where 
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E[w?U] = y 1 {(3y l + 16y 2 +81>g + a(3y l +y 2 +y^) 

+ _L(15y 2 + 48y 2 2 + 243 y 2 + 24y,y 2 + 72y 2 y 3 + 54y 1 y J )) 

SR 2 

2 

E[w? U] = ^ 2 {(y t + 48 y 2 + 81^3) + a(y t + 3y 2 (35 ) 

+ _L_(3y 2 + 240 y 2 + 243y 3 2 + 24y 1 y 2 + 216y 2 y 3 + 18^^} 

SR 2 

2 

E[w?U] = — y 3 ^(yi + 16y 2 + 243 yj + a(y t +y 2 + 3y.,) 

+ J_(3y 2 + 48y 2 + 1215y 3 2 + 8y 2 y 2 + 216y 2 y 3 + 54 y 1 y 3 )) 

SR 2 

Substitution of Eq. (35) into Eq.(34) results in a set of algebraic nonlinear equations for y v y 2 , 
and y 3 . For different excitation levels characterized by a* the foundation modulus a and various 
radii of gyration R, Eq.(34) can be solved numerically. In our study, the standard Levenberg- 
Marquardt algorithm is implemented. 

With y l =E[w i 2 ] obtained, we arrive at the mean square response of the beam as 


3 3 


„ . . mnx . nnx 

E[w 2 (x, 0] - E E sin— sm. 

n«l m*l 


(36) 


L L 

However, the modal amplitudes w m are uncorrelated. Therefore, while the membrane stress 
causes the modal amplitudes to become statistically dependent, it still leaves them uncorrelated. 
Eq.(36) then reduces to 
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E[w 2 (x,t)] = £ y m sin 


2 mux 


m* 1 


(37) 


3. COMPARISON WITH OTHER METHODS 
(1) Fokker-Planck equation method. 

From the solution of Fokker-Planck equation, the exact expression of the probability 
density function is obtained as 


c 


where c is the normalization factor, i.e. 
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at * 
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(38) 
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(39) 




Eq.(38) coincides with that of Herbert (1964) except that we have introduced an additional term 
associated with the elastic foundation. Hence, the modal mean square responses are obtained by 


integration 


? 1 T® r oo 2 I 

E[w m \ = - j. m dw i J_ M '» ex P | 




E " <w - E 

0/1 


N N 
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dw 
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«*i 


r*l 


n* 1 


(40) 
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Eq.(40) must be evaluated by numerical integration. 


(2) Conventional Stochastic linearization Technique. 

By the conventional equivalent linearization method, we obtain 


= n A cOq +au)o + 


2 2 
o) 0 n 2 

AR 1 


E[w^m 2 w 2 ] 


Erf] 


From the equivalent linear system (20), we have Eq.(31). 

Substituting k^ n) in Eq.(42) into Eq.(31) yields 

a 2 = (n* + a)£[w 2 ] 

AR m-i 

(n = 1, 

If only the first three modes are considered, in view of Eq.(27), we obtain 


(41) 


(42) 


o\ = (1 + a ) yi +-L-(Zyl + 4y t y 2 + 9y t y 3 ) 

AR 

Oq = (16 + a)y 2 + J—(y l y 2 + 12y 2 2 + 9y 2 y 3 ) ( 43 ) 

R 

o\ = (81 + a)y 3 + -l-(y l >’ 3 +27 >3 2 ) 

AR 

For specific values of a 0 2 , a and R, one can evaluate E[w 2 ], E[w 2 2 ] and E[w 2 ] through solving 
the set of nonlinear equations Eq.(43). 
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4. RESULTS AND DISCUSSION 

Numerical computations have been performed for the mean square deflection at the 
midspan of the beam for various values of three parameters o 0 2 , a and R. The results from the 
three methods are presented in Figs. 2 through 5. As pointed out by Seide (1975) and as is 
confirmed in equation (37), accurate values of the mean square deflections can be obtained by 
using a three-term approximation, while the number of terms required for accurate stress values 
may be much larger. In the present study, only the first three modal displacements are 
considered, with the emphasis on demonstrating the effectiveness of the new stochastic 
linearization technique. 

As is seen from the differential equation (18), when formally R tends to infinity the effect 
of the nonlinearity disappears. Therefore, magnitude of HR can be viewed as the parameter 
related with the magnitude of nonlinearity. The effect of the foundation stiffness on mean square 
maximum deflection at the midspan of beam is shown in Fig. 2 for one set of levels of 
excitation and nonlinearity. It can be seen that the new equivalent linearization method yields 
more accurate results than the conventional technique. Also, it is shown that the new method 
usually gives greater values of mean square deflection than the exact solution, while the 
conventional method yields values below the exact one. Furthermore, when the stiffness of the 
foundation ^becomes larger (and consequently the parameter a is larger too), both methods tend 
to produce equally accurate results. This conclusion should have been anticipated in view of the 
fact that the system is "more" linear in this case, due to linearity of the Winkler foundation 
model. However, when a is small, the new stochastic linearization method achieves much better 
estimate of the mean square deflection than the conventional technique. The mean square 
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deflections vs parameter R are shown in Figs. 3 and 4; the new method performs much better 
than the conventional stochastic linearization technique for the relatively high nonlinearity of the 
system, i.e. R is from 0 to 1. The effect of strength of excitation on the accuracy of the two 
methods is shown in Fig. 5, where the vertical axis denotes the percentage error of the mean 
square deflection at the midspan of the beam either between the conventional linearization 
method or the new linearization method, and the exact solution; the new approach achieves more 
accurate results than the conventional technique for all the excitation levels. 

To get additional insight into the performance of the proposed method the other set of 
boundary conditions was investigated. Namely, the beam clamped at both ends under the both 
space-wise and time-wise white noise was considered. Fig. 6 portrays the mean square 
displacement calculated by the proposed stochastic linearization method, conventional stochastic 
linearization method and exact solution. An exact solution follows the derivation given in Eqs. 
(14)-(32) except that instead of the sinusoidal mode shape in Eq. (8), the following mode shape 
is utilized 

■qr(x) = cosh y j x - cosy >. x - a(sinh y. x - siny^ x) (^) 

where 


a. 

i 



if m is odd 


if m is even 


(45) 


and the values of y ] are the consecutive solutions of the transcendental equation 
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cosh y- cos y ; = 1 


(46) 


As is seen from Fig. 6 for the clamped beam too, the proposed method results in the mean square 
response which is much closer to the exact solution, than the classical linearization method. 

In both cases of the beams considered the exact solutions were also obtained. The natural 
question arises: How does the proposed method perform when exact solution is not available? 
To answer this question, the additional loading condition was also investigated. Namely, the load 
q(x,t) was represented as a product r(x)q(t). Whereas q(t) was assumed to be weakly stationary 
Gaussian random process namely white noise, r(x) was taken as a deterministic function. 
Specifically r(x) was taken as a constant, representing space-wise uniformly distributed load. 
This representation is valid for the members of relatively short length when the correlation length 
of the excitation is much greater than the length of the beam. For such a loading condition an 
exact solution is unavailable and instead, Monte Carlo simulations should be conducted to check 
the accuracy of the proposed stochastic linearization. Fig. 7 depicts results of such a comparison 
for the simply supported beam. As is clearly seen, the proposed method again exhibits much 
higher accuracy than the conventional linearization technique. 

To sum up, for different boundary conditions and the loading patterns, the suggested 
method is superior to the classical stochastic linearization technique, especially in the important 

high nonlinearity range of the parameters. 
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appendix 


Uncoupledness of Equivalent Linear System 

Under .he assumption .ha, the modal displacements are Normally distributed, one can show tha, 
the equivalent linear system is uncoupled. Indeed, suppose that the equivalent linear system is 

governed by the vector equation 

x (A-l) 

Mw + C>V + Kw = f 


where M and C are diagonal mass and damping matrices, and K is non-diagonal stiffness matrix, 
and fs are independent white noises. New equivalence criterion requires 


E\ 


[ U(w) - —w T Kwf j. = min 




(A* 2) 


where 


W = {tv, H>J - W »} r 

K = [fc] 

The condition that the derivatives of (A 2) with respect to *, equal zero leads to 

E[ww t Kww t ] = 2E[wU(w)w t ] 

For simplicity, let us consider the two-degree-of-freedom system. Suppose 





(A-3) 


(A-4) 


(A-5) 
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Then 


ww T Kww T = 


\w*k u * 2w?w 2 k n + ”> 2*11 + 2w?w?k n + w i w 2 fc 22 

jw 3 w 2 * u + 2w?w?k u + ^ ^w 2 k n + 2w t w 2 k l2 + w 2 fc 22 


w^C/ w Y w 2 U 


E[wU(w)w T ] = 


IW. 


l w 2 J 7 > v 2 2 U 


(A-6) 


Therefore Eq. (A-4) becomes 


( A 1 W - 2 N 


(A-7) 


where 


£[ Wl 4 ] Etw'w 2 ] 


[AJ - 


^[w/w,] 2E[w>j] £[^^2 ] 

£[ w 2 w 2 ] 2 E [ h ' 1 w 2 ] E[w 2 ] 


(A-8) 


ft } = {*U *12 * 22 }’ 


{BJ = {£[w 2 C7] E[w iW 2 U] E[wlU]\ 
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The assumption of zero-mean Gaussian distribution of the modes w„ implies 


In view of Eq. (21), we also have 


E[w?w 2 ] » 0 
E[w^] = 0 


£[WiW 2 17] = 0 

As a result, A x and B 1 turn out to be 


[AJ = 


ftvvj 4 ] 0 E[w?wf] 


0 2E[w?wZ] 0 


E[w 0 £[w 2 4 ] 


(A*9) 


(A- 10) 


(A' 11) 


{fij = {£[^17] 0 E[wtU]} T 

Substituting [AJ, {Bj} in Eq. (A-ll) into Eq. (A-7), one obtains 

k l2 = 0 (A' 12 ) 

This illustrates the uncoupledness of the equivalent linear system in the two-degree-of-freedom 
setting. Analogous proof holds for N > 2. 
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Effect of nonlinearity (R) on mean square deflection at tne mu 
supported beam ( o„ = 10, a =0) (space-wise white noise loading). 



=10, u 
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ipported beam (oj; = 10, a = 1 ) (space-wise white noise loading). 
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at both ends (o 2 0 = 10, a = 1) (space-wise white noise 



, p/pA=0.1 
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Approximate Solution for Random Vibrations 
of Nonlinearly Damped Systems by Partial 
Stochastic Linearization 


I. Elishakoff and G. Q. Cai 


Center for Applied Stochastics Research and Department of Mechanical Engineering 
Florida Atlantic University, Boca Raton, FL 33431-0991, USA 

ABSTRACT The accuracy of the stochastic linearization method is improved by 
the proposed method of partial stochastic linearization, in which only the nonlinear 
damping force in the original system is replaced by a linear viscous damping, 
while the nonlinear restoring force remains unchanged. The replacement is based 
on the criterion of equal mean work, performed by the nonlinear damping force 
in the original system and its linear counterpart. The resulting nonlinear stochastic 
differential equation is then solved exactly, keeping the equivalent damping 
coefficient as a parameter, which can be determined for a specific system by 
solving a nonlinear algebraic equation. 

INTRODUCTION 

There exist very few nonlinear stochastic problems amenable to exact solutions. Therefore, the 
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investigators are resorting to various approximate techniques. One of the most popular 
techniques is the method of stochastic linearization 1_6 . This method has many drawbacks 
however. For example, the tail probabilities for the Duffing oscillator, predicted by the 
approximation may differ form the exact value by a factor as large as 250 7 , and the first 
excursion probability may be in error by several orders of magnitudes 8 . 

Recently, an alternative technique, called "stochastic nonlinearization" has been 
suggested 9 ' 13 , in which one replaces the original nonlinear stochastic differential equation by 
another "close" nonlinear equation, possessing an exact solution. In some cases, however, this 
method is not amenable to closed form solution. 

This study combines the stochastic linearization and nonlinearization techniques. It is 
specially designed for systems with both nonlinear damping and nonlinear restoring force. 
Instead of the classical stochastic linearization technique, where both nonlinear damping and 
nonlinear restoring force are replaced by their respective linear counterparts, here we use only 
a partial linearization, namely, we linearize only the damping. The equation thus obtained is 
amenable to an exact solution. The equivalent damping parameter is obtained by solving a 
nonlinear algebraic equation either analytically or numerically. For the system with nonlinear 
damping but with linear restoring force, the present method coincides with the usual stochastic 
linearization. For the system with both nonlinear damping and nonlinear restoring force, it 
represents a natural generalization of the stochastic linearization. The proposed procedure 
considerably improves the accuracy of the stochastic linearization and yields simple equations to 
determine the desired probabilistic characteristics. 
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BASIC EQUATIONS 


Consider the following nonlinear equation 

X + f(X,X) * g(X) = W(t) (1) 

Here f(X,X) is associated generally with nonlinear damping and g(X) is associated with 
nonlinear restoring force. The excitation W(t) is assumed to be a Gaussian white noise. We 
replace the nonlinear damping force f(X,X) by an equivalent linear one, 

X + + g(X) = W(t) (2) 

Thus, the key issue is to determine the equivalent damping coefficient (3 e for the substituting 
system (2). The criterion for selecting is that the average energy dissipation remains the 
same 11,12 , namely, 

E[Xf(X,X)] = f \E[X 2 ] (3) 

The left and right hand sides of equation (3) represent the average work per unit time performed 
by the original nonlinear damping force and the equivalent linear damping force, respectively. 
It is noted that the stochastic linearization uses the least mean square criterion, namely, 

E { [ f(X, X) - P e A r ] 2 ) = minimum (4) 

which also leads to equation (3). Although the formally identical criterion (3) is used in both 
stochastic linearization and partial linearization, the ensemble averaging in equation (3) is 
performed with different probability densities for these two methods. In the stochastic 
linearization, the probability density is assumed to be Gaussian, while in the present partial 
linearization method, it is generally not Gaussian. 
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We now calculate the left hand side of equation (3). Equation (1) can be written in the 


following Ito differential equation form 


dX l = X 2 dt 

(5a) 

dX 2 = [-f(X lt X 2 ) -g(XJ]dt+y/2nK dB{t) 

(5b) 


where X Y -X, X 2 -X, K is the spectral density of the white noise excitation W(f), and B(t) is a 
unit Wiener’s process. According to the Ito differential rule 14 , 

~Xl = -2 X 2 [f(X lf X 2 ) +g(XJ] +2nK + 2y/2^KX 2 dB(t) (6) 

at 

Ensemble averaging of equation (6) results in 

4- E i X ?] - -2 EiX 2 [f(X v X 2 ) *g^]) +2 kK (7) 

at 

which reduces to 

E{X 1 f(X I ,X,)} = nK (8) 

for the stationary state. Equation (8) leads to a remarkable conclusion: for any oscillator with 

only an additive white noise excitation, the average work done by the damping force per unit 
time depends only on the spectral density of the excitation, regardless the damping mechanism. 
The same conclusion was reached by Kamopp 15 by using a different proof. Hence, from 
equations (3) and (8) 

= E[xjvar\ m (9) 

E[X 2 ] E[X 2 ] 

We will make use of the fact that the exact stationary probabilistic solution of equation 
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(2) is a classical result of the random vibration; it reads 3,4 


p(x,x) = Cexp { -JjL [1L. + j,g(u)du]) 


( 10 ) 


where C is a constant determined from the normalization condition 


C = 


i-i 


J / exp {"-^[-y + j>g( u ) du ] )dxd * 


(ii) 


Equation (10) shows that the velocity of system (2) is a Gaussian random variable, while the 
displacement is not. The probability density (10) can be considered as an approximate one for 
the response of the original system and can be used in equation (8) to yield 


C j j xf(x,x)exp{-J^[?L + (g(u)du])dxdx = kK 


i 


( 12 ) 


Substituting Eq. (11) into Eq. (12), we obtain 



[xf(x,x) -nK] exp {- 


j iK 



X 





= 0 


(13) 


This is a nonlinear algebraic equation for which can be solved analytically or numerically for 
a given system. Once (3, is determined, it can be substituted into equation (10) for the 
approximate joint probability density. 

If the nonlinear damping force depends only on the velocity, namely, / (X,X) -f (X ) , then 
equation (13) is reduced to 
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00 


( 14 ) 


[xf(x) - jiX] exp( ~^l—)dx = 0 
ZtiK 

Equation (14) can often be solved analytically for |3 f . 

It can be seen that if the damping force depends on both the displacement and velocity, 
then the equivalent damping coefficient calculated from (13) is different from that obtained from 
equivalent linearization, because the displacement is not a Gaussian random variable. When the 
damping force is dependent solely on velocity, both methods yield the same equivalent damping 
coefficient from equation (14) since the velocity is Gaussian. However, the present method is 
simpler than the stochastic linearization method since only an equivalent damping coefficient 
needs to be calculated. 

It is noted that Caughey 10 proposed an approximation scheme to replace the original 
nonlinear damping by an energy-dependent nonlinear damping by using the least mean square 
criterion. But his method was only applied to systems with linear restoring force. Lin 13 extended 
Caughey’s method to the case of nonlinear restoring force. However, the determination of the 
replacement nonlinear damping force usually requires numerical calculations. 

As expected, the partial linearization method yields more accurate results than the 
equivalent linearization since it retains one of the characteristics of the original nonlinear system, 
namely the nonlinear restoring force. However, the method is expected to be less accurate than 
the energy dissipation balancing method in which the equivalent damping force is selected from 
a larger class of linear and nonlinear damping forces, rather than from the sub-class of linear 
damping. Therefore, the selected equivalent damping using the energy dissipation balancing 
method is generally "closer" to the original damping. Nevertheless, we pay a penalty of more 
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numerical computations for obtaining more accurate results, when the approximate probability 
density has to be calculated numerically at every point in the energy dissipation balancing 
method, as shown in the illustrative example. One can visualize that for the preliminary design 
of structures, one may need a sufficiently accurate analytical technique such as the proposed 
partial linearization rather than a fully numerical technique. For final design of the structure once 
its parameters are nominally chosen based on partial linearization, one may use a refined analysis 
of the structure vis the dissipation energy balancing method. 

MOMENTWISE CONSISTENCY OF THE PROPOSED METHOD 

It is of interest to note that, in terms of certain statistical moments of the response, the partial 
linearization method is a consistent approximation procedure. Consider M(X l ,X 2 ) -X^X{ where 
i and j are nonnegative integers. An equation for E[M] can be obtained by using Ito’s differential 
rule for M, and then taking the ensemble averaging, resulting in 

4-E[M] ~ *nKE{^i) (15) 

dt dX l dX 2 aX, 2 

The left hand side of (15) is a time derivative of a nth ( n=i+j ) order moment, while the form of 
the right hand side depends on functions / and g. It contains only the nth and lower order 
moments if both / and g are linear; however, it also contains moments of orders higher than n 
if at least one of the / and g functions is a nonlinear polynomial. In either case, equation (15) 
represents a set of relations among moments. 

When the system response eventually attains the stationary state after a sufficiently long 
exposure to Gaussian white noise excitations, the moment equations of the form (15) are reduced 
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to algebraic equations as follows 




dX, 


dX 


nKE[ d —] = 0 

ax 2 


(16) 


Letting M =X V X 2 , X*, X t X 2 and X 2 , we obtain from (16) the following relations: 


£[XJ = 0 (17a) 

ElfiX^XJ+giXJ] =0 (17b) 

E[X t X 2 ] =0 (17c) 

£[*,] -E(X l [f(X v X 1 )+g<X,)}) ( 17d ) 

ElX 2 \J(X v X 2 ) ^(AT,)]! =■ %K (He) 


Equations (17a)-(17e) are satisfied if the true stationary probability density of the system response 
is used when performing the ensemble averaging. They may or may not be satisfied if an 
approximate probability density is used. 

Yet, the use of the approximate probability density (10) guarantees that equations (17a)- 
(17e) are also satisfied provided that g(Xj) is an odd function, and/(A r 1 ,X) * s an °dd function 
of X 2 and an even function of X lf which is true for most physical systems. The validity of the 
above assertion is obvious for equation (17a)-(17c). Equation (17e) is also satisfied by virtue of 
equation (12). Lastly, equation (17d) can be verified by noting that E[X l f(X l ,X 2 )] =0 and 
E\X 2 g(X x )] =£[AT 2 2 ] upon the following integrations by parts: 
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ft x 2 x ' 

c LL x ' gix ' )exp ^ [ -? + J* g(u)du])dx l dx 2 = 


and 


00 00 

C j J x 2 2 exp 
-00 -00 


<-A[£ 

n K l 2 


"l 

+ J>g(d)du])d. 


x, Jx, 


nK 

"E" 


ILLUSTRATIVE EXAMPLE 

For illustration, consider a system governed by 

X + + aX 3 + yX + bX 3 = W(f) 

The stochastic linearization yields 


8 i 1 

p(x,x) = Cexpf- — L(_x 2 + _fcx 2 )l 
nK 2 2 ‘ 1 


where 


P, 


= £ + 
2 


P "\2 


N 


(f) 


3caiK 


and 


* -I + 

' 2 


N 


( 4) 2 


3nX6 


With partial linearization method, only the same (5, is needed to evaluate 
probability density 


( 18 ) 


(19) 


( 20 ) 


( 21 ) 


( 22 ) 


(23) 

"better" 
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( 24 ) 


p(x,x) = C 1 exp[-_?L(i.A: 2 + J_yjc 2 + J-6x 4 )] 
xK 2 2 4 

In this example, the present method is simpler and more accurate than the equivalent 
linearization. 

Let us contrast now the partial and "full" stochastic linearization 1 ' 6,16,17 methods with the 
more accurate energy dissipation method 11 . The probability density is found to be 


P(x,x ) = C,exp 


f (2X -yx 2 - —bx 4 y* dx 
_ P , a r L 2 


. P > _ a f 
xK xKi 


dX 


(25) 


i 1 

J( 2X -yx 2 - —dx 4 ) 7 dx 


where X is the total energy given by 


X = —x 2 + —yx 2 +-Ldx 4 
2 2 4 


(26) 


and 


N N 


(l) 2 +4± -1 , 

6 6 6 


6^0 


A = 


- < 


2 — 


N 


6=0, y #0 


(27) 


The integrals in equation (25) cannot be obtained in closed form if both 6 and y are nonzero. 

By using the method proposed by Caughey 10 and Lin 13 and choosing two coefficients in 
the replacing nonlinear damping force, the approximate probability density is obtained as 
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where X is shown in (26), and the two coefficients c 01 and c 03 satisfy 


m^c, 


1 


02-01 + ( m 04 + Y m 22 + ^ 6 m 4 2 ) C 03 = Mqz + am « 


04 


(29a) 


(m 04 +ym 22 +lbm 42 )c 0l + (m 06 +y 2 m 42 +i-6 2 m 82 +2 Y"» 24 J bm 44 + &ym 62 }c, 


= p(m 04 +ym 22 +ldm 42 ) +a (m 06 + ym 24 +ldm 44 ) 


(29b) 

where m i . = £[li], Equations (28) and (29) can be solved together for c 01 and c 03 . 

The approximate results (25) and (28) may improve the accuracy by introducing a higher 
order term of X. However, numerical computations have to be performed to obtain the 
approximate results. 

In Figs. 1 and 2, the stationary mean square values of the displacement X for system (20) 
are plotted against the stiffness nonlinearity parameter 6 and the damping nonlinearity parameter 
a, respectively. Results computed from both the "full" linearization and partial linearization are 
shown and compared with the Monte-Carlo simulation results. It is seen that considerably higher 
accuracy can be achieved with the present method as compared to the "full" linearization 
procedure. 


CONCLUSIONS 


For a single-degree-of-freedom nonlinear oscillator subjected to an external Gaussian white noise 



excitation, the partial linearization method not only yields more accurate results than those 
obtained from the full linearization method, but also requires less computation. It is also shown 
that the partial linearization method is a consistent approximation scheme in the sense that the 
obtained approximate probability density satisfies certain exact relationships for certain statistical 
moments of the system response. 
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